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One of the main applications of future quantum computers will be the simulation of quantum mod¬ 
els. While the evolution of a quantum state under a Hamiltonian is straightforward (if sometimes 
expensive), using quantum computers to determine the ground state phase diagram of a quantum 
model and the properties of its phases is more involved. Using the Hubbard model as a prototypical 
example, we here show all the steps necessary to determine its phase diagram and ground state 
properties on a quantum computer. In particular, we discuss strategies for efficiently determining 
and preparing the ground state of the Hubbard model starting from various mean-field states with 
broken symmetry. We present an efficient procedure to prepare arbitrary Slater determinants as 
initial states and present the complete set of quantum circuits needed to evolve from these to the 
ground state of the Hubbard model. We show that, using efficient nesting of the various terms 
each time step in the evolution can be performed with just 0{N) gates and Cl(logN) circuit depth. 

We give explicit circuits to measure arbitrary local observables and static and dynamic correlation 
functions, both in the time and frequency domain. We further present efficient non-destructive ap¬ 
proaches to measurement that avoid the need to re-prepare the ground state after each measurement 
and that quadratically reduce the measurement error. 


I. INTRODUCTION 

Feynman envisioned that quantum computers would 
enable the simulation of quantum systems: an initial 
quantum state can be unitarily evolved with a quantum 
computer using resources that are polynomial in the size 
of the system and the evolution time.^ However, it is not 
clear that this would enable us to determine the answers 
to the questions of greatest interest to condensed mat¬ 
ter physicists. The properties most readily measured in 
experiments do not appear to be simple to determine in 
a quantum simulation, and known algorithms for quan¬ 
tum simulation do not return the quantities of great¬ 
est physical interest. In experiments in a solid, one can 
rarely, if ever, prepare a known simple initial state, nor 
does one know precisely the Hamiltonian with which it 
would evolve. Meanwhile, the ground state energy (or 
the energy of any eigenstate yielded by quantum phase 
estimation^^"*) is not a particularly enlightening quan¬ 
tity and, furthermore, given a quantum state it is not 
clear how to determine its long-wavelength ‘universal’ 
properties. Consider, for instance, the following ques¬ 
tions: given a model Hamiltonian for the electrons in a 
solid, such as the Hubbard model, can we determine if 
its ground state is superconducting? Can we determine 
why it is superconducting? In this paper, we show that it 
is possible, in principle, to answer the first question and, 
if it has a clear answer (in a sense to be discussed), the 
second question as well. Moreover, we show that it is not 
only possible in principle but, in fact, feasible to answer 
such questions with a quantum computer of moderate 
size. 

The first step in the solution of a quantum Hamilto¬ 
nian is to map its Hilbert space to the states of a quan¬ 
tum computer. While the Hilbert spaces of systems of 


bosons or spins are mapped most naturally to the states 
of a general-purpose quantum computer, fermionic sys¬ 
tems can also be simulated by using a Jordan-Wigner 
transformation to represent them as spin systems the 
cost can be reduced to log(7V) using additional qubits^ 
or using appropriate ordering techniques discussed below. 
In order to evolve the system on a general purpose quan¬ 
tum computer, we decompose the time-evolution opera¬ 
tor according to the Trotter-Suzuki decomposition®’”; the 
number of time intervals in such a decomposition is de¬ 
termined by the desired accuracy. The evolution for each 
time interval is expressed in terms of the available gates. 
If the Hamiltonian can be broken into k non-commuting 
terms, then the gates can be broken into k sets; within 
each set, the gates can be applied in parallel, and only 
0{k) time steps are needed for each interval. Time evolu¬ 
tion can be used in conjunction with the quantum phase 
estimation^^^ algorithm to find an approximate energy 
eigenvalue and eigenstate of the Hamiltonian. 

The challenge when applying this approach to elec¬ 
tronic systems is that the Coulomb Hamiltonian in sec¬ 
ond quantized form 

N 1 

H — ^ ^ hpqC^^Cq -t- — ^ ^ hpqj-sC^C^CrCs • (1) 

p,q=l p,q,r,s=l 

has k = 0{N‘^) terms for N orbitals. Here the oper¬ 
ators cj) and Cp create and annihiliate an electron with 
spin cr in spin-orbital p (combining orbital index i and 
spin index cr). The quadratic hpq-tevms arise from the 
kinetic energy of the electrons and the Coulomb poten¬ 
tial due to the nuclei (which are assumed to be classi¬ 
cal in the usual Born-Oppenheimer approximation) and 
the quartic hpqrsAerms encode the Coulomb repulsion. 
In order to estimate the ground state energy of a sys- 
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tern of interacting electrons in a large molecule of N or¬ 
bitals, naive estimates indicated that 0{N^) operations 
are needed*^, due to the large number of non-commuting 
terms in the Hamiltonian, but more recent analyses indi¬ 
cate that ~ iV® operations may be sufficient.'^ ^’^2 while 
this makes the simulation of small classically-intractable 
molecules feasible, the scaling is still too demanding to 
perform electronic structure calculations for more than a 
few hundred orbitals. This makes the brute-force simula¬ 
tion of large molecules and crystalline solids impractical. 

One approach to reduce the scaling complexity for the 
simulation of crystalline solids is to focus on effective 
models that capture only the relevant bands close to 
the Fermi energy, and simplify the 0{N^) terms of the 
Coulomb interaction to just a few dominant terms. One 
of the paradigmatic effective models for strongly corre¬ 
lated fermions, discussed in detail in Sec. II, is the square 
lattice Hubbard model, which is a radical simplification of 
the full electronic structure Hamiltonian (I) to arguably 
the simplest interacting single-band model, described by 
the Hamiltonian 

^^Hub = — 

{id) 

+U E E'.".. (2) 

i i 

where = c] ^Ci^cr is the local spin density and rii = 
J2a- total local density. Of all the orbitals in 

a unit cell, we only keep a single orbital (describing e.g. 
the Cu dx 2 -y 2 in a cuprate superconductor), reduce the 
long-range Coulomb repulsion to just the local repulsion 
U within this orbital, and the hybridizations ty = —hpg 
to nearest neighbors on the lattice. The on-site ener¬ 
gies will be called = ha. Usually one can consider a 
translationally-invariant model and then drop the indices 
in tij and since all sites are equivalent. Here we will 
keep them explicitly since they may not all be the same 
during the adiabatic evolution. 

The computational complexity of the Hubbard model 
(2) is much reduced compared to the full electronic struc¬ 
ture problem (I), since only a single orbital is used per 
unit cell instead of dozens needed for the full problem 
and the number of terms is linear in N instead of scaling 
with the fourth power iV^. On a lattice of 20 x 20 unit 
cells this means a reduction of the total number of terms 
from about ~ 10'^ to ^ 10^, which makes such a sim¬ 
ulation feasible on a quantum computer. Furthermore, 
since the 2D Hubbard Hamiltonian can be expressed as 
a sum of five non-commuting terms (the on-site interac¬ 
tion terms and four sets of hopping terms, one for each 
nearest neighbor to a given site) each time step in the 
Trotter-Suzuki decomposition requires only ^ log N par¬ 
allel circuit depth (in fact, constant, if it were not for the 
Jordan-Wigner strings). To achieve this optimal scaling 
we use a Jordan-Wigner transform combined with opti¬ 
mal term ordering and the “nesting” and Jordan-Wigner 
cancellation technique of Ref. 11 . 


The simplifying features of the Hubbard model that af¬ 
ford such advantageous scaling with N are the restriction 
to nearest-neighbor hopping and on-site interactions. Of 
course, a real solid will have long-ranged Coulomb in¬ 
teractions and longer-ranged hopping terms. However, 
the point of a model such as the Hubbard model is not 
to give a quantitatively accurate description of a solid, 
but rather to capture some essential features of strongly- 
correlated electrons in transition-metal oxides. For these 
purposes, a simplified model that does not include all 
of the complexity of a real solid should be sufficient (al¬ 
though, for the case of the cuprate high-Tc superconduc¬ 
tors, the required simplified model may be a bit more 
complicated than the Hubbard model; we focus here on 
the Hubbard model for illustrative purposes). 

A central question in the study of a model of strongly 
correlated electrons is: what is its phase diagram as a 
function of the parameters in the model? The Hubbard 
model has the following interesting subquestion: is there 
a superconducting phase somewhere in the phase dia¬ 
gram? Once the phase diagram has been determined, we 
wish to characterize the phases in it by computing key 
quantities such as the quasiparticle energy gap and phase 
stiffness (superfiuid density) in a superconducting phase. 
To answer these questions we will need to determine the 
ground state wave function for a range of interaction pa¬ 
rameters and densities and then measure ground state 
correlation functions. 

We take the following approach to solve the Hubbard 
model (or related models) on a quantum computer: 

1. Adiabatically prepare an approximate ground state 
|Jtg) of the Hubbard model starting from different 
initial states. 

2. Perform a quantum phase estimation to project the 
true ground state wave function jd^o) from the ap¬ 
proximate state I'I'o), and measure the ground state 
energy Eq = (^'o|iJ|^'o). 

3. Measure local observables and static and dynamic 
correlation functions of interest 

The general strategy for simulating quantum mod¬ 
els has already been discussed previously, starting from 
Feynman’s original proposal.' Abrams and Lloydwere 
the first to discuss time evolution under the Hubbard 
Hamiltonian, but without giving explicit quantum cir¬ 
cuits; in Ref. 14 the same authors discuss how quantum 
phase estimation can be used to project from a trial state 
into an energy eigenstate. Refs. 5 and 15 provided more 
details of how to simulate quantum lattice models, pro¬ 
posed an algorithm to construct arbitrary Slater determi¬ 
nants and discussed measurements of various observables 
of interest. The idea of adiabatic state preparation for 
interacting fermion systems was proposed in Ref. 5 and 
adiabatic preparation has been considered by a number 
of authors'"^'® 

In this paper we go beyond the earlier work in sev¬ 
eral crucial aspects. Besides presenting full details of the 
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necessary quantum circuits for time evolution (Sec. IV) 
and measurements (Secs. VI and VII) we address crucial 
issues that have not received much attention so far. In 
particular, in Sec. Ill we discuss in detail how to adia- 
batically prepare approximate ground states It^g) of the 
Hubbard model starting from different initial states, cap¬ 
turing various proposed broken symmetries, and use this 
to determine the phase diagram of the Hubbard model. 
Having thereby deduced the correct initial state, we can 
cheaply prepare multiple copies of the ground state of the 
system. 

Due to the benign scaling of circuit depth with logfV, 
the ground state of the Hubbard model on large lattices 
can be prepared in very short time of less than a second, 
even assuming logical gate times in the range as we 
discuss in Sec. V. 

We then present, in Sec. VI, details of how to measure 
various quantities of interest, including both equal time 
and dynamic correlation functions. For the latter, we 
present a new approach, based on importance sampling in 
the frequency domain to complement the usual approach 
of measuring in the time domain. 

Since each measurement only gives limited informa¬ 
tion, the ground state has to be prepared many times, 
which leads to a significant total runtime of a hypothet¬ 
ical quantum computer. In Sec. VH, we present new 
approaches to non-destructive measurements of arbitrary 
observables that substantially reduce the runtime. By us¬ 
ing a relatively cheap quantum phase estimation to recre¬ 
ate the ground state after a measurement, we reduce the 
number of times the ground state needs to be prepared 
from scratch. Furthermore their runtime scales as 0(e“^) 
with the desired accuracy e instead of the 0{e~^) scaling 
of the naive approach to measurements. 

For the aficionado of quantum algorithms, the paper 
contains the following new algorithmic ideas. In Sec. V, 
we introduce a variant of quantum phase estimation that 
gives a two-fold reduction in the number of time steps and 
a four-fold reduction in rotation gates compared to the 
standard approach. Section HI presents an approach to 
speeding up adiabatic state preparation by adiabatically 
changing the Trotter time step. In Sec. VI, we present 
an algorithm for measuring dynamical correlation func¬ 
tions by importance sampling of the spectral function in 
the frequency domain. Finally, Sec. VH presents two 
approaches to non-destructive measurements. 


II. THE HUBBARD MODEL 
A. Derivation of the model 

The square lattice Hubbard model is one of the 
paradigmatic models for strongly interacting fermions. 
Although it was originally introduced as a model for 
itinerant-electron ferromagnetism^®’^^, it is now studied 
primarily as a model for the Mott transition and broken- 
symmetry (and, perhaps, even topological) ordering phe¬ 


nomena in the vicinity of this transition (for recent re¬ 
views, see, for instance. Refs. 21 and 22). In particular, 
it has been hypothesized that the square lattice Hubbard 
model captures the important ingredients of the cuprate 
high temperature superconductors. It is a radical sim¬ 
plification of the full Hamiltonian of a solid. A complete 
description of a solid takes the form: 


^ = E 


A. 

2 to . 






a>b 


ZgZbe^ 
\Ra — Rb| 


( 3 ) 


where = 1,...,^ and a,b = l,...Vions; Ne and 
Aions are, respectively, the number of electrons and ions; 
and Zg is the atomic number of ion a and Mg its mass. In 
certain situations, the physics is dominated by electron- 
electron interactions. The ions form a crystalline lattice 
and the vibrations of this lattice (phonons), though quan¬ 
titatively important, do not, it is hypothesized, play a 
major qualitative role. Then we may, in a first approx¬ 
imation, ignore the dynamics of the ions and, thereby, 
obtain a purely electronic Hamiltonian which, in general, 
will take the form of Eq. (1). Explicitly introducing la¬ 
bels for spin, unit cells and orbitals with a unit cell of a 
periodic solid we write this Hamiltonian as 


^ ~ 'y ^ A h.C.^ + ^ ^ 

i,j\a,h\(7 i,a 

2,a 

+ E + . ■ • (4) 

i,j,k,l;a,b,c,d 

Here, i,j label unit cells and a,b label orbitals within 
a unit cell. The local density in orbital a is rii^g = 
coupling constants have the following 
meanings: is the matrix element for an electron to 

hop from orbital b in unit cell j to orbital a in unit cell 
i (possibly two different orbitals within the same unit 
cell); Cg is the energy of orbital a in isolation; Ug is the 
energy penalty for two electrons to occupy the same or¬ 
bital in the same unit cell. This term can be rewrit¬ 
ten equivalently as g Uani,g,^ni^a^i. The prime on the 
summation in the final line indicates that we have omit¬ 
ted the term with i = j = k = I and a = b = c = d, 
which has been separately included as the Ug term. They 
include, for instance, the terms with i = k ^ j = I, 
a = c ^ b — d, which is the density-density interaction 
V^^j^’ni^gUj^b, and the terms with j = I, b = d, which 

is correlated hopping term gCk,cnj^b- These terms 

will be assumed to be much smaller than Ug and will 
be dropped. 

If the differences between the CqS are the largest energy 
scales in the problem, then there will be, at most, a single 
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partially-filled orbital and we can ignore all of the other 
orbitals when studying the physics of the model at tem¬ 
peratures less than \ea — eb\- These energy differences can 
be on the order of eV and, therefore, such an approxima¬ 
tion will be valid over a very large range of temperatures 
and energies. Similarly, if is a large energy scale, it 
may be possible to reduce the model to one with a single 
orbital per unit cell (which, in this case, will be a linear 
combination of the orbitals in Eq. (4)). If the on-site 
interactions U in this orbital (we drop the subscript a 
since there is only a single orbital under consideration) 
are much larger than the interactions with nearby unit 
cells (nearest-neighbor, next-nearest neighbor, etc.), then 
we can drop the latter. Similarly, if the nearest-neighbor 
hopping is much larger than more distant hopping matrix 
elements, then we can drop the latter and write down the 
simplihed model of Eq. (2). This is the Hubbard model. 

In the case of the copper-oxide superconductors, the 
only orbital that is retained is the Cu d^ 2 _y 2 orbital. 
(However, it has also been argued that a 3-band model 
including O Px and Py orbitals is necessary to capture 
the physics of the cuprates. ^^) Usually one can drop the 
indices in tij, Ui and since all sites are equivalent, but 
we will keep them here since during the adiabatic solution 
they will not all be the same. 

Due to the presumed separation of energy scales in (4) , 
the effective model (2) is much simpler than (4) and has a 
much smaller Hilbert space. As we will discuss in detail in 
this paper, this simplification makes simulation of 20 x 20 
or more unit cells feasible on a quantum computer. 


B. The physics of the Hubbard model 

At half-filling (one electron per unit cell) the Hub¬ 
bard model gives a simple account of Mott insulating 
behavior.For U = 0, the ground state is metallic. 
There is a single band and it is half-filled; the Fermi 
surface is the diamond \kx ^ ky\ = 'nja^ where a is the 
lattice constant. However, for [/ > 0, the ground state is 
insulating. While this is true for all U > 0, the physics 
is different in the small U t and U ^ t limits. For 
U t, the system lowers its energy by developing an¬ 
tiferromagnetic order, as a result of which the unit cell 
doubles in size. There are now two electrons per unit 
cell and the system is effectively a band insulator. The 
antiferromagnetic moment N and the charge gap Acharge 
are both exponentially-small N ^ Acharge t for 

some constant c. For large t/, i.e. U ^ t, on the other 
hand Acharge U. In this limit, it is instructive to ex¬ 
pand about t = 0. Then, the ground state is necessarily 
insulating: there is one electron on each site and an en¬ 
ergy penalty of U to move any electron since some site 
would have to be doubly-occupied. For small but non¬ 
zero t/U, this picture is dressed by small fluctuations in 
which an electron makes virtual hops to a neighboring 
site and then returns. As a result of these fluctuations, 
there is an effective interaction between the spins of the 


electrons on neighboring sites. Since an electron can only 
hop to a neighboring site if its spin forms a singlet with 
the spin of the electron there (due to the Pauli princi¬ 
ple), the energy is lowered by such processes if neighbor¬ 
ing spins are antiferromagnetically-correlated. We can 
derive an effective Hamiltonian at energies much lower 
than U by taking into account these virtual hopping pro¬ 
cesses. It takes the form: 

H = (5) 

ihj) 

where J = The 0{t^/U'^) terms are due to vir¬ 

tual processes in which multiple hops occcur If these 
terms are neglected, this is the antiferromagnetic Heisen¬ 
berg model. Its ground state has an antiferromagnetic 
moment of order 1. The temperature scale for the onset 
of antiferromagnetic correlations is ^ J, which is a much 
lower scale than the charge gap ^ U, unlike in the small- 
17 limit, where both scales are comparable.^® When the 
subleading 0{t^/U'^) terms are kept, the system remains 
insulating, but the spin physics may change considerably. 
The phase diagram of such spin models is the subject of 
considerable current research (see, for instance. Refs. 27 
and 28). 

The situation is even less clear when the system is 
“doped”, i.e. when the density is changed from half¬ 
filling. Then, an effective model can be derived to lowest- 
order in t/U. The model is governed by the t — J Hamil¬ 
tonian: 

^ = -EE T h.c.^ -\-J ES.S j+0{t^/U'^) 

{id) {id) 

( 6 ) 

together with a no-double-occupancy constraint. As a 
result of this constraint, the t — J model has a smaller 
per site Hilbert space than the Hubbard model, which 
makes it a much more attractive target for exact diago- 
nalization. However, the t — J model does not capture 
all of the physics of the Hubbard model (especially in 
the regime in which t/U is small but not so small that 
higher-order \nt/U terms can be neglected). 

Any deviation from half-filling will cause the conduc¬ 
tivity to be non-zero in a completely clean system, in 
which is independent of i in Eq. (2). However, in any 
real material, there will be some disorder due to impuri¬ 
ties, and the system will remain insulating for densities 
very close to half-filling. As the density deviates from 
half-filling, antiferromagnetic order is suppressed and 
eventually disappears. If the Hubbard model captures 
the physics of the cuprate superconductors, then super¬ 
conducting order^® with dx^-y^ pairing symmetry must 
occur for a density of 1 — a; electrons per unit cell, with 
X in the range 0.05 ~ x ~ 0.25. Moreover, on the under¬ 
doped side of the phase diagram, 0.05 ~ x ~ 0.15, other 
symmetry-breaking order parameters, such as stripe or¬ 
der, d-density wave order, or antiferromagnetic order may 
also develop (for recent reviews, see, for instance. Refs. 
21 and 22; see, also. Ref. 30 for a recent calculation find¬ 
ing several competing orders). On the other hand, if the 
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Hubbard model is not superconducting in this doping 
range, then superconductivity in the cuprates must be 
due to effects that are missing in the Hubbard model, 
such as longer-ranged interactions, longer-ranged hop¬ 
ping, interlayer coupling, and phonons. 

Thus, as far as the cuprates are concerned, the prin¬ 
ciple aim of an analysis of the Hubbard model is to de¬ 
termine if it has a dx 2 -y 2 superconducting ground state 
for U/t ^ 8 over the range of dopings 0.05 ~ x ~ 0.25. 
If the answer is in the affirmative, then various proper¬ 
ties of this superconducting ground state can be studied. 
For instance, the gap and its dependence on the doping, 
A(a;), can be compared to experimental measurements. 
If the Hubbard model does not superconduct, then dif¬ 
ferent models must be considered, restoring terms that 
have been dropped in the passage to the Hubbard model. 


III. ADIABATIC PREPARATION OF THE 
GROUND STATE 

We start the adiabatic preparation from the ground 
state of some Hamiltonian whose solution is known. For 
the sake of concreteness, we consider two examples. The 
first, in Sec. HI A, is the quadratic mean-field Hamilto¬ 
nian whose exact ground state is the BCS ground state 
for a d-wave superconductor or the analogous Hamil¬ 
tonian for any other ordered state such as striped'^^, 
antiferromagnetic'^^, d-density wave'^^, etc. states. The 
second, discussed in Sec. HIB is the Hubbard model 
with hopping matrix elements tuned so that the system 
breaks into disconnected 2x2 plaquettes. In either case, 
the ground state is known and there is a gap to all excited 
states. 

An essential point is that these initial states can be 
efficiently prepared on a quantum computer. Our ap¬ 
proach to preparing a Slater determinant, discussed in 
Sec. HI A 2, is both deterministic and has better scaling 
than the algorithm previously proposed in Ref. 5. 

We then adiabatically evolve the Hamiltonian into the 
Hamiltonian of the Hubbard model with an additional 
weak symmetry-breaking term (to give a gap to Gold- 
stone modes, a subtlety that is discussed in Section HI A). 
In Sec. IV we review the quantum circuit used to im¬ 
plement time evolution of the Hubbard model and show 
how it can be done with 0{N) gates and a parallel circuit 
depth of only C>(log A) operations per time step. 

This step is formally similar to adiabatic quantum 
optimization'^®, but with a different viewpoint. If no 
phase transition occurs along this adiabatic evolution, 
then we conclude that we guessed correctly about the 
phase of the system (for this value of parameters in the 
final Hubbard model). However, if a phase transition 
does occur, then the gap will close and we will know that 
we guessed incorrectly and can repeat the procedure with 
a different initial state in a different phase. Through a 
process of elimination, we determine the phase of the 
system for a given set of parameters and, by repeat¬ 


ing this process for different Hubbard model parameters, 
we map out the phase diagram. In adiabatic quantum 
optimzation'®, it is assumed that a phase transition oc¬ 
curs and it is hoped that the minimum energy at the 
transition point scales polynomially in system size rather 
than exponentially, thereby allowing the evolution to oc¬ 
cur in polynomial time. In our protocol for solving the 
Hubbard model, however, any gap closing, either poly¬ 
nomial or exponential - corresponding, respectively, to 
a second- or first-order phase transition®' - is an invita¬ 
tion to repeat the procedure with a different initial state 
and correspondingly different annealing path until no gap 
closing occurs. When the initial state is in the same phase 
as the ground state of the Hubbard model, the adiabatic 
evolution can be done in constant time up to subpoly¬ 
nomial factors (see V C for a more detailed discussion of 
errors in adiabatic evolution; while in most of the paper 
we treat the time required for adiabatic evolution as scal¬ 
ing with inverse gap squared, in that section we provide 
a more precise discussion of diabatic errors (i.e., those 
transitions out of the ground state that are due to a fi¬ 
nite time scale for the evolution) and further show how 
it is possible to improve the scaling with gap by using 
smoother annealing paths. 


A. Adiabatic Evolution from Mean-Field States 

1. Mean-Field Hamiltonians 


Let us suppose that we wish to test whether the ground 
state of the 2D Hubbard model is superconducting for 
some values of t, 17 and density. We can do this by seeing 
if the Hubbard model is adiabatically connected to some 
simple superconducting Hamiltonian. It need not be real¬ 
istic; it merely needs to be a representative Hamiltonian 
for a superconductor. We can take the BCS mean-field 
Hamiltonian: 


rrMF _ 
^DSC — ~ 




(ij) 


- E - AA) + h-c-’ 

{i,j) 


where A“ = A/2 for i = / ± x and A“ = —A/2 
for i = j±y. This is a mean-field dx 2 -y 2 superconductor 
(DSC) with fixed, non-dynamical superconducting gap of 
magnitude A. 

The ground state of this Hamiltonian can be written in 
the form of a Slater determinant by making a staggered 
particle-hole transformation on the down-spins: cl^ —> 
Then the Hamiltonian takes the form: 


£VMF 

^DSC 


EE^b (4. 


cr^jjCr A Cj (jCj 




- E(-i)'^b (4tS7 - +h-c., 

(i.j) 


( 8 ) 
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Under this transformation, the number of up-spin elec¬ 
trons per unit cell is particle-hole transformed: riii —)■ 
1 — riii while the number of up-spin electrons is un¬ 
changed; equivalently, the chemical potential for up-spin 
electrons is flipped in sign, Consequently, 

the number of down-spin electrons is not equal to the 
number of up-spin electrons. In addition, the supercon- 
ducting pair field A“- ^ has become transformed into a 
staggered magnetic field in the Sx direction with a dx 2 -y 2 
form factor (a d-wave spin-density wave, in the terminol¬ 
ogy of Ref. 33). Thus, the ground state of this quadratic 
Hamiltonian is a Slater determinant of the form: 

1^) = + ^kcl,+Q^^)|0) (9) 

k 


Here, the momenta k range over the Brillouin zone, 
— - < kx V < where a is the lattice constant, and 
the functions Uk, Uk are given by: 



Cfc 


y/el + 


( 10 ) 


where Ck = —2t{cos kxQ + cos kyu) and Afc = A(coskxa — 
coskya). In Sec. HI A 2, we show how to prepare this 
ground state. For now, we will simply assume that this 
ground state can be prepared efficiently. In computing 
the properties of this state, it is important to keep in 
mind that physical operators have been transformed ac¬ 
cording to cl^ —>■ (—l)*Cj^ and to use the appropriate 
transformed operators. 

We now adiabatically deform this Hamiltonian into a 
weakly-perturbed Hubbard model. There are two rea¬ 
sons why our final Hamiltonian is a weakly-perturbed 
Hubbard model, rather than the Hubbard model itself: 
(1) in the absence of long-ranged Coulomb interactions, 
(which are neglected in the Hubbard model) a super¬ 
conducting ground state will have gapless Goldstone ex¬ 
citations, and (2) a dx 2 -y 2 superconductor has gapless 
fermionic excitations at the nodes of the superconduct¬ 
ing gap. Although these gapless excitations are impor¬ 
tant for the phenomenology of superconductors (e.g. the 
temperature dependence of the superfluid density and the 
structure of vortices), they are a nuisance in the present 
context. Hence, we weakly perturb the Hubbard model in 
order to give energy gaps to Goldstone modes and nodal 
fermions. However, this can be done weakly so that we 
do not bias the tendency towards superconductivity (or 
lack thereof). If the Hubbard model superconducts, then 
the scale of the dx 2 -y 2 superconducting gap should be 
~ xJ, where (n^) = I — x is the average occupancy of a 
site and J = AP/U is the superexchange interaction. (On 
the other hand, if the Hubbard model does not supercon¬ 
duct, then the superconducting gap in the cuprates must 
be determined by some other energy scale.) Hence, we 
apply a fixed dx 2 -y 2 superconducting gap that is much 


smaller than xJ in order to pin the phase of any su¬ 
perconducting order that may develop. In addition, we 
apply a small imaginary dxy superconducting gap to give 
a small gap to the nodal fermions. Neither gap scales 
with system size. 

More concretely, we adiabatically evolve the Hamilto¬ 
nian according to: 


Hoscis) = (1 - s)-f^DSC + 




(bi> 


E -N? ( 


- AA 


^ -I- h.c. (II) 


Here, w is a small dimensionless parameter that deter¬ 
mines the magnitude of a small pair field that remains 
throughout the adiabatic evolution. At s = 0, it simply 
increases the superconducting gap by a factor of 1 -I- w. 
At s = I, it applies a small pair field oc w to the Hub¬ 
bard model. If the Hubbard model does not have a su¬ 
perconducting ground state, then w ^ 0 will induce a 
small superconducting gap oc w. However, if the Hub¬ 
bard model does have a superconducting ground state, 
then w will pin the phase of the superconducting gap, 
which would otherwise fluctuate. The magnitude of the 
gap will be essentially independent of w for small w. If 
the Hubbard model is superconducting, then the system 
will remain superconducting all the way from s = 0 to 
s = 1 and the gap will remain open. As long as s is 
varied slowly compared to the minimum gap - which, 
crucially, is independent of system size - the system will 
remain in the ground state. The final line of Eq. (II) is 
a small imaginary second-neighbor superconducting gap. 
Here, AJ^ = uA/2 for j = k± {x + y), AJ^ = uA/2 for 
j = A: ± (x — y), and it is a small dimensionless number 
that is independent of system-size. Such a term opens a 
small gap oc u at the nodes (assuming that the Hubbard 
model does not have a superconducting ground state that 
spontaneously breaks time-reversal symmetry). We take 
uA xJ so that this term opens a gap at the nodes 
but is otherwise too small to affect the development of 
superconductivity away from the nodes. 

In summary, Eq. (II) is a path through parameter 
space from a simple superconducting Hamiltonian with 
Slater determinant ground state to a Hubbard model that 
has been perturbed by small terms that give a gap to or¬ 
der parameter fluctuations and nodal excitations but oth¬ 
erwise do not affect the ground state, as may be checked 
by decreasing the magnitude of these perturbations. 

A similar strategy can be used to test whether the Hub¬ 
bard model has any other ordered ground state, such as 
antiferromagnetic order. In this case, our starting Hamil- 
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tonian is, instead: 


Slater determinant state is a state 


rrMF 

^AF 






( 12 ) 


We then evolve along a path analogous to Eq. (11), 
but with the symmetry-breaking terms in the last two 
lines replaced by a weak staggered magnetic field. If the 
Hubbard model has an antiferromagnetic ground state 
(for some values of t, U, x) then the gap will not close 
along this path. However, if the actual ground state of 
the Hubbard model is superconducting for these values 
of the couplings, then the gap will become 
where 0 is the dynamical critical exponent of the transi¬ 
tion (i.e. the gap will close, up to finite-size effects) at 
the point along the evolution at which superconductivity 
develops. There will be a non-zero staggered magnetiza¬ 
tion throughout since we never fully remove the staggered 
magnetic field, but there will be a sharp signature at the 
point at which superconductivity develops’^"*. The pres¬ 
ence of a staggered field, which is kept small, may slightly 
shift the onset of superconductivity but will have no other 
effect. 

This strategy can allow us, in some of the circum¬ 
stances in which it has a clear answer, to address the 
question: why does the 2D Hubbard model have a su¬ 
perconducting ground state? In particular, it is possible 
that the superconducting ground state of the Hubbard 
model necessarily has some other type of order coexisting 
with superconducting order, especially in the underdoped 
regime. These secondary orders play a role in some theo¬ 
ries of the cuprate superconductors. If such an order were 
present in the Hubbard model, evolution from a purely 
superconducting initial state would necessarily encounter 
a phase transition (at the point of onset of this additional 
order). Therefore, we can determine not only if the Hub¬ 
bard model has a superconducting ground state, but also 
whether the ground state exhibits a secondary type of 
long-ranged order over some range of dopings. We will 
return to the question “why does the Hubbard model 
superconduct” when we discuss correlation functions. 

The virtue of using these mean field Hamiltonians 
as starting points for adiabatic evolution is that their 
physics is fully understood. In particular, their ground 
states are Slater determinants, which we can efficiently 
prepare, as we explain in the next section. 


2. Preparing Slater Determinants 

In this section, we explain an efficient, simple quantum 
algorithm to prepare Slater determinants. This is then 
used to prepare the Hubbard ground state by adiabatic 
evolution from [/ = 0 to C/ 7 ^ 0 . 

The standard algorithm is in Ref. 5. This algorithm 
suffers from some drawbacks. Given a projector p, a 


'i'sD(p) = nf^,&]|0), (13) 

where | 0 ) is the no particle vacuum state, and each op¬ 
erator b\ is a linear combination of a], where the a] are 
fermion creation operators on a given site i. We write 
= Tri=i4PiP where P 4 are complex scalars and are 
the entries of P, n is the number of orbitals, and Ne is 
the number of electrons. Given any matrix p, we can find 
a matrix P such that 


p = PP^ (14) 

Then, using this matrix P in the definition of the state 
gives the state '^sd{p)', note that this definition fixes 
'^Sd{p) up to a complex phase. The algorithm in Ref. 5 
proceeds by noting that Um = exp(i|(&m -|- 6 j„)), acting 
on the vacuum, produces the state f)|„| 0 ), up to a phase. 
Therefore, if the are suitably orthogonalized, succes¬ 
sively acting with the unitary Um produces the desired 
state. However, no explicit method for implementing this 
unitary Um is given; if the unitary is implemented by a 
Trotter method, this will be extremely inefhcient. We 
explain a simpler technique. 

Let us first assume that p is real; the extension to 
the complex Hermitian case will be straightforward. The 
idea is first to prepare the Slater determinant state 

^SD(po) = nf^ia]|0). (15) 

This corresponds to the case that p is an n-by-n diagonal 
matrix, with ones in the first N^. entries and zeroes else¬ 
where. This state has one electron on each of the first 
fVe states. It is a product state that can be prepared in 
linear time (indeed, simply initialize each of the first 
spins to up and the other spins to down). Then, we act 
on this state ^ sd{p) with a sequence of Givens rotations. 
A Givens rotation, to borrow the language of matrix di- 
agonalization, is a rotation matrix that acts on two rows 
or columns at a time. For any given pair, i, j, define the 
operator Rij{0) by 

Rij{0) = exp( 0 c|cj- — h.c.). (16) 

This is a simple two-qubit gate (up to Jordan-Wigner 
strings) and is, therefore, more straightforwardly imple¬ 
mented than the unitary Um in the standard method 
described above. Define the n-by-n orthogonal ma¬ 
trix ri j{9) by its matrix elements as follows. Let 
{rij{9))k,k = I for k ^ i,j and let (rij{B))kg = 0 if 
k ^ I and k ^ i,j or ii k ^ I and I yf i,j. Then, 
let {rij{9))i^i = {rij{9))jj = cos( 6 >) while {rij{9))ij = 
= sin( 6 »). Then, 

R^A0)'^SD{p) = ^SD{r^A0)pr^.Ji0)~"), (17) 

where again the equality holds up to the phase ambiguity. 

Now, given some desired target 4 > 5 £)(p), we claim that 
one can always find a sequence of at most nNg Givens 


rotations that will turn the matrix p into the matrix po- 
Thus, by inverting this sequence of Givens rotations and 
applying that sequence to the state 'I'sd(po)j we succeed 
in constructing the desired state. The Givens rotations 
will be ordered such that the Jordan-Wigner strings can 
be cancelled as in Ref. 11, so the algorithm takes time 
0{Nf^n). To show our claim that the sequence of Givens 
rotations exists, note that the effect of a given rotation 
which transforms p to rij{9)prij(0)~^ can be achived by 
transforming P to rij{9)P. Consider a left singular vec¬ 
tor of P with singular value equal to 1; we can find a 
sequence of at most n Givens rotations that rotate that 
vector to be equal to the vector ( 1 , 0,...). As a result, all 
the other singular vectors now have amplitude only on 
sites 2,...,n. Find another sequence to transform some 
other left singular vector with singular value equal to 
1 into the vector (0,1,0,...). Continue until all singu¬ 
lar vectors with nonzero singular values have amplitudes 
only on the first Ne sites. 

The extension of this procedure to the complex case is 
simple. Rather than having the Givens rotation Rij{9) 
depends only a single angle 9, we need to perform a rota¬ 
tion exp{zclcj — h.c.) for some complex number z. All the 
time estimates remain unchanged up to constant factors. 

However, for many systems a much shorter sequence of 
Givens rotations can be found. Using the same idea as in 
the fast Fourier transform, if n is a power of 2, there exists 
a sequence of n log 2 (n) Givens rotations to transform an 
n-by-n matrix to the momentum basis. Thus, we can 
produce ground states of a system of free fermions on a 
periodic lattice in time O{n^\og{n)). This can be easily 
extended to handle the case of a periodic system with a 
unit cell larger than a single site. Assume there are states 
|i, x) where i labels an atom in a unit cell and x labels a 
lattice vector. We transform to a momentum basis \i, fc), 
and then for each k we initialize some Slater determinant 
in the unit cell in momentum space. If the unit cells have 
size 0 ( 1 ), the time to initialize each cell is 0(1) so the 
total time is still 0{n? log(n)). 

A further possible optimization is that if W > 'n/2, we 
can instead use Givens rotations to bring the holes to the 
desired states, rather than the electrons, taking a time 
0({n — Ne)n) rather than 0{N^n). As an interesting 
application of this, consider creating the ground state 
of A^e = 3 electrons in n = 4 sites, using a hopping 
Hamiltonian with the sites arranged on a ring. This can 
be done by initializing a state with electrons on the first 
3 sites and then apply a total of three different Rij(9) 
operators. 

3. Numerical Results 

We applied the Slater determinant preparation pro¬ 
cedure to the Hubbard model on 8 sites arranged in two 
rows of 4 sites. We started at [/ = 0 and Ci = 0. All hori¬ 
zontal couplings were set equal to 1 , while to avoid having 
an exact ground state degeneracy the vertical couplings 


were doubled in strength to 2 (note that this requirement 
depends on the filling factor; for example, at half-filling 
we would instead prefer to remain at vertical coupling 
equal to 1 to avoid degeneracies). The ground state of 
this model is a Slater determinant of free fermion single¬ 
particle wavefunctions that we could prepare with a total 
of 14 Givens rotations for both spin up and spin down. 
This was done by initializing particles in three out of four 
sites in the top row of sites, and holes in the remaining 
sites. Then, 3 Givens rotations were used to create a 
uniform superposition of the hole in various states in the 
top row. Finally, 4 Givens rotations were used to create 
a uniform superposition between top and bottom rows. 
This was done for both spins, giving 2 x 7 = 14 rotations. 

We then annealed to t/ 7 ^ 0 while reducing the verti¬ 
cal coupling to unit strength. The results are shown in 
Fig. 1. As seen, increasing the annealing time increases 
the success probability, until it converges to 1. The an¬ 
nealing times required to achieve substantial (say, 90%) 
overlap with the ground state are not significantly differ¬ 
ent from the preparation approach discussed next based 
on “joining” plaquettes. 

Note that for this particular model, we know the 
ground state energy with very high accuracy from run¬ 
ning an exact diagonalization routine so we know whether 
or not we succeed in projecting onto it after phase esti¬ 
mation. The system size is small enough that calculating 
ground state energies is easy to accomplish on a classical 
computer using Lanczos methods. Of course, on larger 
systems requiring a quantum computer to simulate, we 
will not have access to the exact energy. In this case, it 
will be necessary to do multiple simulations with differ¬ 
ent mean-field starting states as described at the start 
of this section; then, if the annealing is sufficiently slow 
and sufficiently many samples are taken, the minimum 
over these simulations will give a good estimate of the 
ground state energy. Once we have determined the op¬ 
timum mean-field starting state and annealing path, we 
can then use this path to re-create the ground state to 
determine correlation functions. 


B. Preparing of d-wave resonating valence bond 
states from plaquettes 

1. d-wave resonating valence bond states 

An alternative approach to adiabatic state prepara¬ 
tion has previously been proposed in the context of ul¬ 
tracold quantum gases in optical lattices. In this ap¬ 
proach, d-wave resonating valence bond (RVB) states 
are prepared. RVB states are spin liquid states of a 
Mott insulator, in which the spins on nearby sites are 
paired into short-range singlets. RVB states were first 
introduced as proposed wave functions for spin liquid 
ground states of Mott insulators'^^, and then conjec¬ 
tured to describe the insulating parent compounds of the 
cuprate superconductors.^® In more modern language. 
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FIG. 2. Labeling of the sites in a plaquettes 
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FIG. 1. Success probability as a function of total anneal¬ 
ing time for an annealing path in the Hubbard model that 
begins at 17 = 0 and vertical hopping equal to 2 and ends 
at 17 = 2.5t (blue) or 17 = 4t (red) and all hoppings equal 
to 1. The initial state is a Slater determinant state (of free 
fermion single-particle states). Annealing is linear along the 
path. The a;-axis represents the total time for the anneal. The 
y-axis represents the probability that the phase estimation 
routine returns an energy that is within 10“® of the ground 
state energy. This probability is obtained by sampling runs 
of the phase estimation circuit, giving rise to some noise in 
the curves. 

they are described as Z 2 topologically-ordered spin liq¬ 
uid states.The idea was that, upon doping, the sin¬ 
glet pairs would begin to move and form a condensate of 
Cooper pairs. Although the insulating parent compounds 
of the cuprates are antiferromagnetically-ordered, it is 
nevertheless possible that the cuprate superconductors 
are close to an RVB spin liquid ground state and are 
best understood as doped spin liquids. The RVB sce¬ 
nario has been confirmed for t-J and Hubbard models 
of coupled plaquettes’^®’^® and ladders (consisting of two 
coupled chains)'*^’"^^ and remains a promising candidates 
for the ground state of the Hubbard model and high- 
temperature superconductors."*^ 

Instead of starting from mean-field Hamiltonians, we 
build up the ground state wave function from small spa¬ 
tial motifs, in this case four-site plaquettes, which are 
the smallest lattice units on which electrons pair in the 
Hubbard model.'*® Following the same approach as orig¬ 
inally proposed for ultracold quantum gases, we prepare 
four-site plaquettes in their ground state, each filled with 
either two or four electrons.*^ These plaquettes then get 
coupled, either straight to a two-dimensional square lat¬ 
tice, or first to quasi-one dimensional ladders, which are 
subsequently coupled to form a square lattice. 

2. Preparing the ground states of four-site plaquettes 

The first step is preparing the ground state of the Hub¬ 
bard model on 4-sites plaquettes filled with either two or 
four electrons. We start from very simple product states 


O-O.o-o 

0 2 4 6 

FIG. 3. Goupling of two plaquettes 


and adiabatically evolve them into the ground states of 
the plaquettes. 

To prepare the ground states of plaquettes with 
four electrons we start from a simple product state 
Cq ^|0), using the labeling shown in Fig. 2. 

Our initial Hamiltonian, of which this state is the ground 
state has no hopping {Uj = 0 ), but already includes the 
Hubbard repulsion on all sites. To make the state with 
doubly occupied state the ground state we add large on¬ 
site potentials £2 = €3 = 2U -b 3t on the empty sites. To 
prepare the ground state of the plaquette we then 

1 . start with = 0, Ui = U and a large potential 
on the empty sites (here £2 = €3 = e) so that the 
initial state is the ground state, 

2 . ramp up the hoppings tij from 0 to t during a time 

Ti, 

3. ramp the potentials Ci down to 0 during a time T 2 . 

Time scales Ti = T 2 ~ 10t“* and e = U -\- At are suf¬ 
ficient to achieve high fidelity. As we discuss below, 
it can be better to prepare the state quickly and then 
project into the ground state through a quantum phase 
estimation than to aim for an extremely high fidelity in 
the adiabatic state preparation. To prepare a plaque¬ 
tte with two electrons we start from the product state 
Cq^cJ^| 0) and use the same schedule, except that we 
choose the non-zero potentials on three sites of the pla¬ 
quette (ei = £2 = £3 = e)- 
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FIG. 4. Probability of preparing ground state as a function of 
times Ti, T 2 ,T 3 in units of , following the annealing sched¬ 
ule described above. During time T3, a uniform potential e 
was ramped from f to 0 on the plaquette with two electrons 
to break the symmetry between plaquettes. 


S. Coupling of plaquettes 

After preparing an array of decoupled plaquettes, some 
of them with four electrons and some with two electrons 
depending on the desired doping, we adiabatically turn 
on the coupling between the plaquettes. As a test case 
we use that of Ref. 17 and couple two plaquettes (see 
Fig. 3), one of them prepared with two electrons and 
one with four electrons. 

Naively coupling plaquettes by adiabatically turning 
on the intra-plaquette couplings t 24 and fas will not give 
the desired result, since the Hamiltonian is reflection 
symmetric but an initial state with two electrons on one 
plaquette and four on the other breaks this symmetry. 
No matter how slowly the anneal is done, the probability 
of preparing the ground state will not converge to 1. We 
thus need to explicitly break the reflection symmetry by 
either first ramping up a small chemical potential on a 
subset of the sites or - more easily - not completely ramp¬ 
ing down the non-zero values when preparing plaque¬ 
tte ground states. Consistent with Ref. 17 we find that a 
time T 3 « is sufficient to prepare the ground state 

with high fidelity. 

See Fig. 4 for success probabilities depending on times 
Ti to ramp up hopping in a plaquette, T 2 to ramp down e 
in a plaquette, and T 3 to join the plaquettes. The specific 
couplings in these figures are an initial potential 8t on the 
empty sites in the plaquette with four electrons, chemical 
potential 8t on the empty sites in the plaquette with two 
electrons and chemical potential t on the occupied site in 
the plaquette with two electrons. See Fig. 5 for a plot of 
the spectral gaps for an example annealing schedule. 

Going from two plaquettes to the full square lattice 
is straightforward: we bias some plaquettes with small 
values of to break any reflection and rotation symme¬ 


15 



FIG. 5. Energy levels during annealing, y-axis is energy, x 
axis is time. From time 0 to time 16 (in units of t~^) the 
hopping is ramping up, time 16 to time 36 the non-uniform e 
in each plaquette is ramping down (leaving a uniform e = 1 
on the plaquette with two electrons), time 36 to time 76 is 
joining plaquettes. All levels remain non-degenerate except 
for a degeneracy that appears at the end of the anneal. After 
time 76 plaquettes are separated, as described below. 


tries, and then switch on inter-plaquette hoppings, and 
switch off the bias potentials. The time scale here is a- 
priori unknown. According to Ref. 17 the ground state 
on ladders can be prepared within a relatively short time 
T « 2001“^ — 500t“^. In two dimensions we expect that 
similar time scales may be sufficient unless we encounter 
a quantum phase transition to a different phase. 


4- Decoupling plaquettes 

As proposed in the context of analog quantum simula¬ 
tion, the pairing of holes on plaquettes can be tested 
by adiabatically decoupling plaquettes. We start by 
preparing one plaquette with four electrons and a second 
plaquette with two electrons. We then couple them as 
described above to end up with a system with two holes 
relative to the half hlled Mott insulator. After adiabati¬ 
cally decoupling the two plaquettes again we measure the 
number of electrons on each plaquette. If holes bind in 
the ground state, they will end up on the same plaquette, 
and we measure an even electron number per plaquette. 
If they don’t bind each hole will prefer to be on a differ¬ 
ent plaquette to maximize its kinetic energy and we will 
end up with three electrons per plaquette. 

As shown in Fig. 6 , the probability of finding three 
electrons on a plaquette goes to zero for U = 2t and 
U = At when increasing the time Tg over which we ramp 
down the inter-plaquette hopping. This indicates pairing 
and is consistent with the known critical value of 17 « 4.5t 
for pairing on four-site plaquettes. On the other hand for 
U = 6t and U = 8t we see an increase in the probability 
of finding three electrons per plaquette, indicating the 
absence of pairing. 
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FIG. 6. Shown is the probability P33 of finding three electrons 
on each of the plaquette after decoupling two plaquettes with 
a total of six electrons as a function of the time Ts (in units 
of t~^) over which the plaquettes were nearly adiabatically 
separated. A decrease of P33 indicates pairing of holes in the 
ground state of a plaquette. 


C. Trotter Errors and Accelerating adiabatic 
preparation by using larger time steps 

In practice, it suffices for our annealing preparation 
to obtain a state with sufficiently large overlap onto the 
ground state that phase estimation will then project onto 
the ground state with a reasonably large probability. We 
can thus tolerate some errors in our annealing path. For 
this reason, we annealed using a larger time step, before 
doing the phase estimation with a smaller time step. 

In fact, this approach of annealing at a larger time 
step can be turned into an approach which uses a larger 
time step but still produces a state with very large over¬ 
lap with the ground state. Note that the larger Trotter- 
Suzuki step means that we in fact evolve under a modified 
Hamiltonian, proportional to the logarithm of the unitary 
describing a single time step. As the time step decreases, 
this Hamiltonian becomes closer to the desired Hamil¬ 
tonian. It is then possible to evolve with a larger time 
step, following the adiabatic path in parameter space, 
and then at the end of the path one can “anneal in the 
time step”, gradually reducing the time step to zero. 

In our simulation of the anneal using joining two pla¬ 
quettes, the system was rather insensitive to Trotter er¬ 
rors. We used a relatively large time step, equal to 
0.25t“^ for the annealing, and used a much smaller time 
step for the phase estimation that made the errors negli¬ 
gible there. 

A time step of 0.05t“^ for phase estimation yielded 
relative errors below 10“^ and 0.003t“^ yielded relative 
errors below I0“®. The larger time step used for the 
adiabatic preparation prevents us from obtaining 100% 
overlap with the ground state even in the limits of very 


Name 

Symbol 

Matrix Representation 

Hadamard 



I/V2 l/y/2 
I/V2 


Y basis change 

s 


1 -1 

i-H 

t-H 

1 _1 


Phase gate 


im 



1 0 

0 


Z rotation 


RM 



giS/2 Q 

0 


controlled not (CNOT) 



10 0 0 

0 10 0 

0 0 0 1 

0 0 10 



TABLE I. An overview of the quantum gates used in this 
paper. Note that Y gates do not represent Pauli Y ; rather, 
they represent a Clifford operator that interchange Y and Z 
spin orientations. The phase gate T{6) can be replaces by an 
Rz{0) rotation, up to a global phase that can be kept track 
of classically. 


long annealing times; however, as seen we still obtain 
very high overlaps in the range of 80% — 96%. 


IV. IMPLEMENTING UNITARY TIME 
EVOLUTION OF THE HUBBARD MODEL 

A. Circuits for the individual terms 

For the Hubbard model, we need to implement uni¬ 
tary evolution under two types of terms, hopping terms 
and repulsion terms cl^^Cp^crcl^^Cq^c,- These 
terms are a subset of those needed for earlier elec¬ 
tronic structure calculations*^. In addition we will need 
to implement unitary evolution under a pairing term 

To establish notation we summarize In Table I the 
gates used in the quantum circuits shown in this paper. 


1. The repulsion and chemical potential terms 

The repulsion term is an example of what is called an 
Hpqqp term and evolution under this term is given in Fig. 
7. A chemical potential term, cj,g.Cp_cr, which we call an 
Hpp term, is not needed for the final Hamiltonian but is 
useful for annealing. The circuit for this term is shown 
in Fig. 8. Since all the repulsion terms in the Hubbard 
model and all the chemical potential terms commute with 
each other, they can be applied in parallel, thus needing 
only 0{N) gates and 0(1) parallel circuit depth. 
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FIG. 7. Quantum circuit to implement time evolution 
for a time step 6 under the term Hpqqp = rip^crriq^^i = 
'-p.o-Cq ^iCq,a'Cp,a- Top and bottom lines represent qubits for 
p, a and q, a'. 


im 


FIG. 8. Quantum circuit to implement time evolution for a 
time step 6 under the chemical potential term = cj, o-Cp^a. 

2. The hopping and pairing terms 

The unitary evolution under the hopping term, that we 
call an Hpq term, is given by the circuit in Fig. 9, and has 
also been considered previously. It also will be useful to 
be able to implement evolution under a super-conducting 
pairing term, Acj^ +h.c.. This term is similar to the 
hopping term, since both are quadratic in the fermionic 
operators. The hopping and pairing term circuits will 
be very similar to each other, with the changes involv¬ 
ing changing some of the basis change gates. The pair¬ 
ing term will be needed for the adiabatic evolution from 
the BCS mean-field Hamiltonian suggested in subsection 
III A 1. Also, it will be needed to measure superconduct¬ 
ing correlations present in the Hubbard state: in this case 
the measurement will be of the correlation function be¬ 
tween two different superconducting pairing terms on two 
pairs of sites separated from each other. Measurement 
will be discussed more later, but our general procedure 
for measurement will require knowledge of the unitary 
which implements evolution under the term. 

For the case A real, we wish to implement unitary evo¬ 
lution by exp[—-I- h.c.)]. Let us write Xp ,j to 
denote the Pauli X operator on the qubit corresponding 
to spin-orbital p, a, and similarly Yp^a- and Zp^a- denote 
Pauli Y and Z operators on that spin-orbital. Using the 
identity that c'^ = (1/2)(A + iY) and c = [1/2){X — lY), 
up to phase factor associated with the Jordan-Wigner 
strings, we wish to implement 

Q 

axp[ ? —(Ap fj Yp fjYq^tyt'^Z (f^) 

where Zjw denotes the product of the Pauli Z operators 
on the spin-orbitals intervening between p,u and q,(j'. 
In general, we wish to implement this unitary operation 
controlled by some ancilla. In fact, this circuit is the same 


as the circuit used to implement the term cj, -bh.c., 
up to a sign change on the second spin rotation, as shown 
in Fig. 10 . 

For A imaginary, we wish to implement 

Q 

(i^P[-i^iXp,aYq^o- + yp,aXq^qr')Zjw]- (19) 

This circuit is similar to the other pairing and hopping 
circuits, except that in the first half of the circuit, we 
apply a Hadamard to the first qubit and a Y basis change 
to the second, and in the second half we apply the Y 
basis change to the first qubit and the Hadamard to the 
second. 


B. Optimizing the cost of the Jordan Wigner 
transformation 

In order to implement the unitary time evolution in the 
Hubbard model, we use a Jordan-Wigner transformation 
to turn the model into a spin model. The depth of the 
circuits described in the previous subsection will depend 
upon the particular ordering of spin-orbitals chosen for 
this Jordan-Wigner transform. We choose to order so 
that all qubits corresponding to spin up appear first in 
some order, and then all qubits corresponding to spin 
down. Since the hopping terms preserve spin, this sim¬ 
plifies the Jordan-Wigner strings needed. For a given 
spin, we order the sites in the “snake” pattern shown in 
Fig. 11. 

With this snake pattern, we group the hopping terms 
into four different sets, as shown. The terms in any given 
set all commute with each other. The two horizontal sets 
of terms have no Jordan-Wigner string required, since 
they involve hopping between neighboring sites with the 
given ordering pattern, and so these terms can all be 
executed in parallel. The vertical terms do not commute 
with each other. However, they nest as in Ref. (11) so 
that all vertical terms in a given set can be executed with 
0{N) gates using a depth 0{Vn), where we assume that 
the geometry is roughly square so the the number of sites 
in a given row is Vn. 

For periodic boundary conditions, we must also handle 
the terms going from left to right boundary or from top 
to bottom boundary. The terms going from left to right 
boundary can be executed in parallel with each other. 
Each term naively requires a depth proportional to the 
length of a given horizontal row to calculate the Jordan- 
Wigner string. Assuming a square geometry, this length 
is 0{'/N). We can use a tree to calculate the fermionic 
parity of the sites in the middle of the string, reducing 
the depth to log( A) if desired. For the terms from top to 
bottom edge, there are Jordan-Wigner strings of length 
0{N), but these terms nest and so can all be executed in 
parallel, reducing the depth to (!I(log(A)). The fermionic 
parity on the sites on the middle rows (those other than 
the top and bottom row) can be calculated using a tree, 
again in depth 0{\og{N)). The total depth can then be 
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FIG. 9. Quantum circuit to implement time evolution for a time step 6 under the hopping term c], o.Cq,cr + cl^^Cp^a- Top and 
bottom lines represent qubits for p, a and q, a. 
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FIG. 10. Quantum circuit to implement time evolution for a time step 9 under the pairing term + Cq,a'Cp^a- Top and 

bottom lines represent qubits for p, a and q, a'. 



FIG. 11. Snake pattern for ordering of sites. Arrow shows 
order of sites. Dashed and solid lines represent hopping terms. 
The hopping terms are grouped into four different sets, two 
sets of vertical and two sets of horizontal hoppings as shown 
by different colorations and dash patterns. 

(!I(log(A^)-\/]V). In fact, using a tree to compute parities, 
the total depth can be reduced to polylogarithmic in iV. 

Another way to reduce the cost of the Jordan-Wigner 
transform in this geometry is in Ref. 43. This method 
does require additional qubits. 


V. FINDING THE GROUND STATE BY 
QUANTUM PHASE ESTIMATION 

The annealing algorithm generates a state with good 
overlap with the ground state if the rate is sufficiently 
slow. V C discusses optimized paths that increase the 


overlap. However, in fact it is not necessary to obtain 
a perfect overlap with the ground state; if the overlap 
with the ground state is reasonably large, then we can 
apply quantum phase estimation to the state and have 
a reasonably large chance of projecting onto the ground 
state. Since the time for the simplest annealing path 
scales inversely with the gap squared, while the time to 
use quantum phase estimation to measure energies more 
accurately than the gap (and thus to determine whether 
or not one is in the ground state) scales inversely with 
the gap, up to logarithmic corrections, using quantum 
phase estimation to project an approximate state onto 
the ground state may be preferred. The gap will depend 
upon the annealing path; since we use a symmetry break¬ 
ing field until nearly the end, the gap only closes at the 
end of the anneal. In this section we first give a simple 
improvement to quantum phase estimation that leads to 
a reduction in depth. We then give some numerical esti¬ 
mates for annealing times for free fermion systems (the 
simplest case where one can estimate annealing times for 
a large system numerically). 


A. Reducing the number of rotations in quantum 
phase estimation 

While this general setting of using adiabatic prepara¬ 
tion and quantum phase estimation has been considered 
many times before, we briefly note one simplification of 
the quantum phase estimation procedure. If the gate 
set available consists of one- or two-qubit Clifford oper¬ 
ations, then this simplification provides a slightly arger 
than twofold reduction in depth; more importantly, it 
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provides a fourfold reduction in the number of arbitrary 
rotations, which are expected to be the most costly op¬ 
erations to implement. This simplification is not in any 
way specific to the Hubbard model, and could be ap¬ 
plied in other settings, such as in quantum chemistry. 
The usual approach to quantum phase estimation is to 
build controlled unitaries to implement the unitary evo¬ 
lution (by using the ancilla to control the unitaries de¬ 
scribed above). In this case, if the ancilla qubit is in 
the |0) state, then no quantum evolution is implemented, 
while if the ancilla qubit is in the |1) state, then quan¬ 
tum evolution hy U = exp{—iHt) is implemented (up to 
Trotter-Suzuki error). Then, quantum phase estimation 
estimates the eigenvalues of this unitary U. The sim¬ 
ple modification that we propose is that instead, if the 
ancilla qubit is in the |0) state, then we implement evo¬ 
lution by exp(iHt/2), while if the ancilla qubit is in the 
|1) state, then we implement evolution by exp(—iHt/2). 
Since the difference between these two evolutions is the 
same exp{—iHt), then the same phase estimation proce¬ 
dure on the ancilla will return the identical result for the 
energy (again up to Trotter-Suzuki error). 

This approach means that we have then only to imple¬ 
ment unitary evolution for half the time (i.e., t/2 rather 
than t), so assuming the same time step for each Trot¬ 
ter step means that the number of Trotter steps required 
is halved. However, to analyze whether or not this im¬ 
proves the gate depth, we need to determine how the gate 
depth might change in a given Trotter step. First, let us 
review the usual method of making circuits controlled to 
implement phase estimation, and then we will explain 
how to do it for the modification considered here. As an 
example, let us consider a circuit to implement a hopping 
term, as in Fig. 9. The usual technique (see references 
before) to make this circuit controlled (so that if the an¬ 
cilla is |0) then the identity operation is performed while 
if the ancilla is |1) then evolution under the hopping term 
is performed) is to modify the two Rz{d) gates, making 
them both controlled by an ancilla. A controlled rota¬ 
tion by angle 9 implements no rotation if the ancilla is in 
the |0) state and implements a rotation by angle 9 if the 
ancilla is in the |1) state. Different quantum computing 
architectures may make different elementary gates avail¬ 
able. If we have access to arbitrary Z rotations by angle 
9, and to Clifford operations, then the controlled rotation 
by angle 9 is implemented by rotating the target by an¬ 
gle 9/2, then applying a controlled NOT from the ancilla 
to target, then rotating the target again by angle —9/2, 
then again applying a controlled NOT from the ancilla 
to target, so that the net rotation is 9/2 — 9/2 = 0 if 
the ancilla is in the |0) state and is 0/2 -|- 0/2 = 0 if the 
ancilla is in the |1) state. Thus, the controlled rotation is 
implemented by doing two arbitrary angle rotations, and 
additionally applying two Clifford operations, for a total 
of four gates given the gate set mentioned above. 

A crucial point is that the only gates which need to 
be made controlled are the Rz{9) gates. The other 
gates in Fig. 9 do not need to be modified in any way. 


The property still holds when we consider our proposed 
method of implementing evolution by either exp(iHt/2) 
or exp{—iHt/2), depending on the ancilla. In this 
method, we need to modify the Rz (0) gates to implement 
a rotation by either a positive or negative angle; we term 
this a “uniformly controlled rotation”: a rotation by an¬ 
gle +9 if the ancilla is in the |0) state and a rotation by 
angle —9 if the ancilla is in the |1) state. This can be 
done by applying a controlled NOT from the ancilla to 
target, then rotating the target, then again applying a 
controlled NOT from the ancilla to target. This requires 
one arbitary angle rotation and two Clifford operations. 

Thus, we find that the total number of gates re¬ 
quired to implement the appropriately controlled version 
of Fig. 9 has been reduced by 1, as one fewer arbitrary 
rotation is required, while the same Clifford operations 
are required. So, the depth of a single Trotter step is 
slightly reduced. Since the number of Trotter steps has 
been halved, we find indeed a slightly larger than twofold 
reduction in depth, as claimed. Further, the number of 
arbitrary rotations is reduced by four. We have described 
this contruction for the hopping term; however, for any 
term, the usual technique is to replace arbitrary single 
qubits rotations by controlled rotations while we suggest 
instead replacing them with uniformyl controlled rota¬ 
tions. 


B. Annealing Time for Free Fermions 

To make some estimate of the cost of adiabatically 
preparing a good approximation to the ground state, it 
is worth considering the case of free fermions. For any 
system of free fermions with translation symmetry and 
k sites per unit cell, the annealing problem decouples 
into many different problems, one for each Fourier mode, 
each problem describing annealing of a fc-dimensional 
subspace. For k = 2, this is a Landau-Zener problem. 
For a single such Landau-Zener problem, the annealing 
time is expected to scale as 1/A^ where A is the gap of 
the Hamiltonian. However, we have many Landau-Zener 
problems in parallel, and we want all Fourier modes to 
reach the ground state. The problem is least serious if 
the Fermi surface consists just of isolated points (such 
as in one dimension or for a semimetal in higher dimen¬ 
sions) as there are only a few Fourier modes with small 
gap. For systems with a Fermi surface, the density of 
states at the Fermi surface is larger, meaning that there 
are more modes with small gap. However, the probabil¬ 
ity of a diabatic transition for a Landau-Zener system is 
exponentially small in the ratio of the gap squared to the 
Landau-Zener velocity. Hence, even if there are a large 
number of modes with gap A, this should require only a 
logarithmic slowdown in the adiabatic velocity to have a 
small probability of a diabatic transition. However, one 
must also worry about effects due to the endpoints of the 
evolution. See, however, V C for a detailed discussion of 
improved annealing strategies that overcome some of the 
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effects of the endpoints. 

To put some numbers onto these generalities, we con¬ 
sidered a system of spinless free fermions at half-filling in 
one dimension with periodic boundary conditions (this 
of course is the most advantageous case with only two 
modes with minimum gap). For each Fourier mode, 
we compute the evolution using numerically exact tech¬ 
niques (this is a two-dimensional quantum mechanics 
problem); to determine the probability of the entire sys¬ 
tem being in the ground state, we multiply the probabil¬ 
ities for each mode. We follow the annealing path 

Hs= (^cla+i + h.c)j + s ^ {c\ci+i + h.c.y (20) 

i even i odd 


starting at a fully dimerized state at s = 0 and moving 
to a uniform state at s = 1, where (ci) creates (an¬ 
nihilates) a spinless fermion on site i . We followed a 
uniform annealing schedule and used a time discretiza¬ 
tion with negligible Trotter error. We chose sizes which 
were not a multiple of 4 to avoid having an exact zero 
mode. 

For a system 1002 sites, with an annealing time t = 
1000 , the probability of ending in the ground state is 
0.54.... Thus, the expected time to end in the ground 
state is roughly 2000 plus an roughly additional 2 phase 
estimation steps on average. Reducing the annealing 
time to t = 500, the probability of ending in the ground 
state is 0.306..., taking slightly less expected anneal¬ 
ing time, but more phase estimation steps on average. 
Increasing to 2002 sites, the situation is worse. For 
t = 1000, the probability of ending in the ground state 
is 0.13..., so the expected time is over 7500, plus over 
7 phase estimation steps on average. Because the sys¬ 
tem is translationally invariant with period 2 for all s, 
the probability of ending in the ground state is a prod¬ 
uct over momentum vectors of the probability of ending 
in the ground state in each momentum vector. Each of 
these probabilities can be calculated by evolving a 2-by- 
2 matrix. In all cases, the dominant contribution to the 
smallness of the probability to end in the ground state 
was from the lowest energy mode. Thus, small changes 
in the low energy density of states can lead to dramatic 
changes in the success probability. 

For 1002 sites, the lowest eigenvalue of the free fermion 
hopping matrix is 0.0188..., so the gap A of the Hamil¬ 
tonian is twice that. Hence, in this case, the phase esti¬ 
mation time could be much smaller than the annealing 
time. 

Since rather large Trotter time steps of order 0.25 are 
sufficient to obtain reasonable overlap with the ground 
state, we find that about 10^ steps are sufficient to pre¬ 
pare the ground state. Using parallel circuit each step 
can be performed with a small circuit depth that in¬ 
creases only logarithmically with the system size. Hence 
a parallel circuit depth of about one million gates is suf¬ 
ficient to prepare the ground state. 


C. Improved Annealing Paths 


There are broadly two strategies for reducing the er¬ 
ror in the annealing path. The typical objective in such 
optimizations is to hnd a function, /(s), that satisfies 
/ : [0,1] !-)■ [0,1] and /(O) = 0 and /(I) = 1 such that 
the diabatic leakage out of the groundstate caused by 
Udiab ■■= Eo{s)ds^^-^T H{f{s}}ds jg minimized for 

a fixed value of T. Here we use the convention that 
T is the time ordering operator, T is the total evo¬ 
lution time allotted and s = t/T is the dimensionless 
time, and the phase factor is included to ensure that 
limr-i-oo Udiab = Uadiab is well defined. 

The first strategy is known as local adiabatic 
interpolation.'^'^ The idea of local adiabatic interpolation 
is to choose /(s) to increase rapidly when the eigenvalue 
gap of H{f{s)) is large and increase slowly when the gap 
is small. This strategy can substantially reduce the time 
required to achieve a hxed error tolerance, but does not 
improve the scaling of the diabatic errors in the state 
preparation, which scale as 0{1/T) for Hamiltonians that 
are sufficiently smooth. 

The second strategy, known as boundary cancella¬ 
tion is often diametrically opposed to local adiabatic 
optimization."*® The idea behind it is instead of mov¬ 
ing at a speed designed to minimize diabatic transitions, 
you choose a path designed to maximize cancellations 
between them. There are several strategies for exploit¬ 
ing this, but the simplest strategy to to choose /(s) such 
that dg f{s)\s=o.i = 0 for A: = 1,..., m. Then if /(s) is 
at least to -I- 2 times differentiable and H{f{s)) does not 
explicitly depend on T then the error in the adiabatic 
approximation is at most*® 


2 


max 

S 


\\dT+^Hif{s))\\ 

A(s)m-l-2ym-|-l 


-hO(l/T™+2) ^ 


( 21 ) 


where A(s) is the gap at the given value of s. Equa¬ 
tion (21) therefore implies that if 


^ foy^(l-yrdy 


( 22 ) 


then the upper bound on the asymptotic scaling is im¬ 
proved from 0(1/T) to 0(l/r'"+*). Eurthermore, if a 
diabatic error of S is desired then it suffices to take 


T = 0 


( max 


\\dT+^H{fism^A ^ 
A(s)*+^(5^ ) 


(23) 


This implies that even in the limit of small 6 the cost of 
adiabatic state preparation is not necessarily prohibitive 
because to can be increased as S shrinks to achieve im¬ 
proved scaling with the error tolerance. Eurthermore, the 
scaling of the error with the minimum gap becomes near- 
quadratically smaller than would be otherwise expected 
(for 6 sufficiently small). 

The above argument only discusses the scaling of the 
evolution time under the assumption that to is fixed. If 








16 


m grows as S shrinks then there may be m-dependent 
prefactors that are ignored in the above analysis. If the 
Hamiltonian is analytic in a strip about [0,1] and it re¬ 
mains gapped throughout the evolution then Lidar et al. 
show that, for fixed T, the error grows at most polyno- 
mially with m.'^® Calculus then shows that the optimal 
value of m scales logarithmically with the error tolerance 
and hence such contributions are subpolynomial. 

These two strategies are often mutually exclusive when 
the avoided crossing does not occur near the boundary 
because if H{s) is analytic then H{s) will have to be 
larger away from the boundary to compensate for the fact 
that it is nearly constant in the vicinity of the boundary. 
This can cause the Hamiltonian to move rapidly through 
the avoided crossing, which increases the annealing time 
required. However, the minimum gap often occurs at the 
end of the evolution when transforming from the initial 
Hamiltonian (such as the free Fermion model, plaque- 
tte preparation, or other tractable approximation) to the 
interacting theory. This means that the two strategies 
are often compatible here and that boundary cancellation 
methods will often be well suited for high-accuracy state 
preparation. Further adiabatic optimizations could also 
be achieved by altering (22) to approximate the local- 
adiabatic path in the region away from from s = 0,1. 

We are of course concerned about the evolution time 
required to perform Uadiab ~ Udiab since it dictates the 
resources needed to prepare the ground state ipo within a 
fixed error by adiabatic evolution starting from the free 
fermion ground state V’q/. A related problem is the prob¬ 
lem of implementing the projector onto the ground state, 
Pq = IV'o)(V'ol within a fixed error; we will wish to be able 
to do this for a method described in VH B. Given the pro¬ 
jector we have Pq = Uad^abPo^Ula^^^,; 

we describe later how to measure Pq and so by conju¬ 
gating this measurement by Uadiab we can measure Pq. 

At first glance, the use of a Trotter decomposition may 
appear to be problematic for adiabaticity because the 
time-dependent Hamiltonian that describes the decom¬ 
position is discontinuous. Such discontinuities are not 
actually problematic, as can be seen by using the trian¬ 
gle inequality 

mad^ab-UZ,)\^^{,f)\ (24) 

^ \\Uadiab ~ C^Ziafcll 
< \\Uadiab - UdiabW + \\Udiab - 

where is a quantum circuit approximation to the 

finite time evolution C/diab- This shows that the error in 
the state preparation is at most the sum of the adiabatic 
error and the error in simulating the adiabatic evolution. 
We have a similar bound for error in measuring Pq: 

\\>-iadiab ^0 ^adiab ^ diab^O ^ diabW 

< 2 {\\Uad^ab - Ud^ab\\ + \\Ud^ab " G|fb||) • (25) 

In order to make the entire simulation error at most d, 
it suffices to set both sources of error to 5/2. This cre¬ 
ates an issue because the cost of simulation scales at least 


linearly with T, which in turn scales with \/5. There 
are a wide array of different Trotter-based formulas that 
can be used in the simulation but in all such cases the 
Trotter number (and in turn the simulation cost) scales'*' 
for integer fc > 0 as which using (21) is 

/ ~(2k+m. + 2) \ 

O I 5 2 '‘("'+i) J For any fixed m, the value of k that leads 

to the best scaling with 5 can found by numerically opti¬ 
mization. However, if we approximate this optimal value 
by taking 2k = m for even m then the cost of the Trot- 
terized simulation is 0(5“^/’”), which leads to subpoly¬ 
nomial (but not polylogarithmic) scaling of the total cost 
of simulation in 1/6 if 2k = m and the Hamiltonian is 
sufficiently smooth because the cost of a Trotterized sim¬ 
ulation grows exponentially with k unlike the error in the 
adiabatic approximation. 

This shows that the cost of implementing Pq or prepar¬ 
ing 'ijjQ within error 5, Tq, is sub-polynomial in 1/5. Since 
C>(e“* log(l/e)) repetitions of this circuit are needed 
and since errors are subadditive it follows that taking 
6 cx e^/log(l/e) suffices to make the total error at most 
e. Thus the cost of implementing Pq scales as 
where the o(l) term is the amalgamated costs of the adi¬ 
abatic preparation and the logarithmic term from phase 
estimation. 

It also should be noted that high-order Trotter-Suzuki 
formulas may not be needed in practice as numerical ev¬ 
idence suggests that small Trotter numbers typically suf¬ 
fice for the small Hubbard models that we have consid¬ 
ered. This means that low-order methods may often be 
preferable to their higher-order brethren. Tight error 
bounds for the second order Trotter-Suzuki formula are 
given in the appendix. 


VI. MEASURING OBSERVABLES 

In this section we will discuss the measurement of in¬ 
teresting physical observables. The total energy, which 
could be measured by phase estimation, is the least inter¬ 
esting quantity. We will instead focus on densities and 
density correlations, kinetic energies. Green’s functions 
and pair correlation functions. 

Such computations will allow us to understand the 
physics of the ground state. For instance, if the ground 
state is superconducting, then we can compute the corre¬ 
lation function of the pair field by the methods described 
in this section. Its long-distance behavior determines the 
condensate fraction. We can also compute the expecta¬ 
tion value of the kinetic energy; its variation with doping 
can help determine whether the system becomes super¬ 
conducting as a result of kinetic energy gain compared 
to the non-superconducting state. Finally, more com¬ 
plex correlation functions, which would be very difficult 
- if not impossible to determine in experiments - could 
determine how these superconducting properties vary in 
response to perturbations that enhance or suppress ef¬ 
fects such as spin or charge fluctuations. 
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A. Local observables and equal-time correlations 

1. Measuring the density, double occupancy and spin and 
density correlations 

The densities are trivially measured by measuring 
the value of the qubit corresponding to the spin-orbital 
{i, a). Similarly double occupancies can be deter¬ 

mined by measuring two qubits, while density correlation 
functions 

n^rij = -I- (26) 

and spin correlation functions 

StSj = - rij-;) (27) 

can be determined directly by measuring the values of 
four qubits. 

By simultaneously measuring all qubits in the compu¬ 
tational basis (eigenvectors of the Pauli-Z operator) we 
can determine density, double occupancy and all spin and 
density correlations. As an example, suppose we have 
recorded this sequence of meaurement outcomes over sev¬ 
eral runs (for each run, for every qubit we measure the 
Pauli-Z operator and we record the measurement out¬ 
comes). Given this data, suppose we now wish to esti¬ 
mate a correlation function such as 

for a pair of sites i ^ j. This is equal to the expectation 
value of 

+ 1 -\- 1 
2 2 ■ 

Given the outcomes of the measurements of 
we simply add 1 to each measurement, multiply the re¬ 
sults, and divide by four. We then average this over runs. 
For example, if we consider a two-site Hubbard model at 
half-filling and U = 0, we will find that the two measure¬ 
ment outcomes Zi^^ and Zj^^ are perfectly anti-correlated 
(since there is only one electron with spin up which has 
equal probability to be on either of the two sites) and so 
the average of the product will not equal the product of 
the averages. 

2. General strategy for other observables 

Next, we note that there is a general strategy used 
to measure the expectation value of any unitary oper¬ 
ator U, assuming that we can build a circuit that im¬ 
plements a controlled version of this unitary, controlled 
by some ancilla. Namely, we apply a one-bit phase es¬ 
timation using the phase estimation circuit of Fig. 12. 
This is a standard trick; see for example Fig. 9 in Ref. 
15. Since we have circuits that implement unitary evo¬ 
lution under various Hamiltonian terms, this enables us 
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FIG. 12. General phase estimation circnit to compute the 
expectation value of any nnitary which can be given as a 
controlled black box. The Z{6) produces a rotation abont 
the Z axis by angle 6-, by varying this, real and imaginary 
parts of the expectation value can be measured. 


to meaure these terms. For example, to measure a term 
c\,iyCq,a: we measure the expectation value of the unitary 
exp[—i9{cjj ,„.Cq^cr+h.c.)]. Since the operator cjj ,„.Cq^a+h.c. 
has eigenvalues —1,0,+1, the most efficient results are 
obtained from the phase estimation when we choose 
9 = TT (which perfectly distinguishes in a single mea¬ 
surement between eigenvalue 0 and eigenvalues ±1) or 
9 = 1 ^ j 2 (which perfectly distinguishes the case of eigen¬ 
value -1-1 from the case of eigenvalue —1). 

However, for the observables we consider, there in fact 
are much simpler ways of measuring the correlation func¬ 
tions. We give two different strategies. The first strat¬ 
egy involves replacing rotations in some of our unitary 
gates with measurements and we call it the “stabilizer 
strategy”; the second introduces a new gate called an 
“FSWAP”. 


3. The Stabilizer Strategy 

The stabilizer strategy is a method for measuring ob¬ 
servables of the form exp[—i0(Oi-|-O2 + ---)] or of the form 
O 1 +O 2 + ■■■■, where each Oi is a product of some number 
of Pauli operators, and [Oi,Oj] = 0. This form includes 
many operators of interest to us, including terms in the 
Hamiltonian such as the kinetic energy as well as other 
terms such as pair correlation. We call this the stabilizer 
strategy because of our use of products of Paulis which 
commute with each other; no assumption is made that 
any state is a stabilizer state. 

If we can measure each Oi , then we succeed in measur¬ 
ing the desired operator. Since they commute, we may 
measure them in any order without needing to recreate 
the state after measurement. As each Oi is a product of 
Paulis, there is some unitary in the Glifford group which 
maps it onto a Pauli Z operator on some given qubit. 
Hence, we can measure that Oi by applying that Glifford 
unitary, then doing a Z-basis measurement, and finally 
undoing the Glifford unitary. Applying this procedure to 
terms in the Hamiltonian, for which we have previously 
given circuits to implement exp[— i0(C)i -I-O 2 + ■•■)] using 
sequences of controlled Z-basis rotations conjugated by 
Clifford gates, the measurement circuit amounts to re¬ 
placing each controlled Z-basis rotation in the evolution 
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circuit with a Z-basis measurement. For example, if we 
apply this procedure to the Hpq unitary evolution shown 
in Fig. 9, we arrive at the circuit shown in Fig. 13. Sum¬ 
ming the value of the two Z-measurements (treating the 
answers as ±1 for each measurement) and then divid¬ 
ing by two (this division occurs because the controlled 
unitaries rotate by an angle equal to half the coupling 
strength) gives a measurement of the expectation oper¬ 
ator of the hopping operator. In fact, we find that the 
same measurements are needed to measure the pairing 
operator for A real, since as we have noted, the pairing 
operator for A real is implemented by the same unitary 
as the hopping operator, up to a sign change in the sec¬ 
ond controlled rotation. Thus, the two measurements we 
have given simultaneously measure the two commuting 
operators + h.c. and c|, + h.c. (to measure 

the pairing operator, one must instead consider the dif¬ 
ference of the two Z measurements). 


giving the unitary: 

/ 1/V2 0 0 1/^2 \ 

-1/V2 0 0 1/V2 , , 

0 1/V2 1/V2 0 ■ ^ ^ 

\ 0 - 1/^2 1 / V 2 0 

An example circuit containing several measurements 
in parallel is shown in Fig. 14. 

We may then measure the Z eigenvalues of the two 
qubits and call them si and S 2 - If we measure si = — 1, 
there is one particle on the two sites. We then have as 
estimators for the Green’s function cjj^^Cq^a' + „,Cp,a- 

S2Ss,,-i (32) 

and for the kinetic energy —tpq (cjj ,^Cq^a- + the 

estimator 


4- The FSWAP Strategy for the kinetic energy and Green’s 
functions 

For the kinetic energies we will have to measure one- 
body Green’s functions cl^Cj^cr + To do this we 

first swap qubit j with its neighboring qubits until it 
is neighboring qubit i in the normal ordering and is at 
one of the positions k = i ± 1. Taking into account the 
fermionic nature of the electrons we need to use a new 
fermionic swap gate FSWAP with matrix 


tpqS2Ssi,— l (^^) 

Note that if we have N qubits we can do N/2 measure¬ 
ments simultaneously, as long as no qubit is involved in 
two measurements. 

Conversely, if the first qubit is -1-1, there are either 0 
or 2 particles on the bond. In the presence of a pair¬ 
ing term, which breaks particle number conservation we 
can measure non-trivial values for the pair field operator 


/ 1 0 0 0 \ 

0 0 10 
0 10 0 
V 0 0 0 -1 / 

The circuit expressed in standard gates is: 


(28) 


X 

T 

X 


X- 


X 




(29) 

We next change the basis to be able to better mea¬ 
sure the kinetic energy on the bond which is then 
~tpq (cp a^q,<r + a^p,cr) ■ This basis change is a simple 
circuit: 


5. The FSWAP Strategy for the pair correlation functions 

Pair correlation functions in their most general 
form are built from the general 2-body terms 
cl c,Cj -a'^k,-a'Cpa- This expectation value can be mea¬ 
sured after fermion-swapping the qubits to be adjacent 
and then performing two Green’s function measurements. 
We first use FSWAP gates to move the qubits i, k, j, 
and I to four adjacent positions which we will label 1 
to 4. We then apply the same circuit as for the Green’s 
function measurements to both pairs (1, 2) and (3,4) and 
measure the Z component of each qubit m (=1,2,3, or 
4) as Sm- The expectation value is then expressed as 

{c\^a^k,aclj^^iCl^cr') = (■S 2 S 4 ^si11^83, —l) 

Care must be taken if two indices are the same, e.g. if 
i = k in above term. In that case things simplify since 

(34) 

and the measurement just becomes 



(35) 

Similarly if * = A: and j = I we get 
(30) 


( 36 ) 
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FIG. 13. Hpq Measurement 



FIG. 14. Example of a circuit for parallel measurements of 
a number of Green’s function values between different sizes 

and the measurement is straightforward. 

For the Hubbard model we are interested in singlet 
pairing and the specific terms we want are: 

^ ~ - Cfc.;Ci,^). (37) 

Multiplying this out and sorting by spin we get 

Xj,ki = “2 

, (38) 

which can be measured by individually measuring all 
terms. 

Note that even though there may be different 
choices for the four indices, we do not need to measure 
all of these four-point functions. As pairs are expected 
to be tightly bound in the Hubbard model, choosing i 
close to j and k close to I, but choosing the pairs {i,j) 
and (fc, 1) at as large a distance as possible on a given lat¬ 
tice is sufficient to check for long-range pair correlations. 
In fact, we can measure -I- h.c (or, — h.c) 

for N/2 pairs (i,j) simultaneously because the operators 


commute; then, averaging the result over pairs extracts 
the long wavelength pair correlation. To measure these 
with minimum depth, one can use nesting strategies sim¬ 
ilar to those discussed for measuing hopping terms in 
Sec. IV B. 


B. Dynamic correlation functions and gaps 

1. Dynamic response in the time domain 

The measurement of dynamic correlation functions 

CA,B{t) = {iPo\A{t)B{0)\il;o) (39) 

for time-independent Hamiltonians such as the Hubbard 
model are typically presented not in the time domain but 
after a Fourier transform in the frequency domain: 

/ OO 

dtCAAty (40) 

-OO 

The standard way of measuring them has been explained 
in several references, most explicitly in Ref. 15. For 
unitary operators A, B, the product is a 

unitary, and can be measured using the circuit of Fig. 12. 
One simple improvement is to note that it is not neces¬ 
sary to control the evolutions and it suffices 

to control A\B] then, if the ancilla bit is |0) so that the 
operators A\B are not done, the product jg 

equal to the identity as desired. The uncontrolled time 
evolution will require only half as many arbitrary angle 
rotations as the controlled time evolution for many gate 
sets (see Sec. V A for discussion of implementing con¬ 
trolled time evolution). 

A further improvement is to note that the final evo¬ 
lution just gives an overall phase acting on the 

ground state, and so can be replaced by Z rotation of 
the ancilla without being implemented, giving another 
factor of two reduction in depth. A final improvement 
combines the idea in Ref. 15 of controlling one opera¬ 
tion and anti-controlling another with one of the ideas 
in Sec. V A of implementing either forward or backward 
evolution in time. We prepare the ancilla in the state |-|-). 
If the ancilla is |0), we implement and if the an¬ 
cilla is |1) we implement The measurement of 
































































































20 


the ancilla in the X basis gives the desired expectation 
value. The controlled time evolution, by either or 

^-itH /2 depending on whether the ancilla is |0) or |1), 
can be done as described in Sec. V A. This requires half 
the number of time steps as that required implement the 
uncontrolled evolution by assuming that the same 
timestep is used for both evolutions; the only additional 
cost is an additional CNOT gate for each arbitrary angle 
Z rotation in the evolution and since these are Clifford 
gates the cost is much smaller than the cost of an arbi¬ 
trary angle rotation for many gate sets. 


2. Dynamic response in the frequency domain by phase 
estimation 


A dynamic correlation function of an operator with its 
adjoint, such as 

CA{t) = {tjjo\A\t)A{0)\ipo) (41) 

can be simplified after a Fourier transform into the fre¬ 
quency domain: 


Sa{oj) 



dte^^*CA(t) (42) 

n 

n 


^|(V'„|A|^o)P<5(cc-(A„-Ao)). 


Instead of performing a “simple sampling” and mea¬ 
suring CA(t) for all times t and then Fourier-transforming 
it and hoping to get these delta-functions resolved a new 
and better approach is to do “importance sampling” and 
measure the energies of the eigenstates \n) with eigenen- 
ergies En directly with the weight |(^/>„|A|'!/)o)P with 
which they appear in above sum. 

This can be achieved by phase estimation. If we ap¬ 
ply phase estimation to the state A\'\j)Q) then eigenstate 
|'0n) will be picked just with the weight A | (■!/;„ |A|'(/)o)P 
and a histogram of the measured energies thus directly 
measures the dynamic structure factor SA(to). 

The retarded correlation function: 


XA(t)=0(i)(^o|[At(t),A(O)]|V^o) (43) 

can then be obtained as follows. First, we write the 
retarded correlation function in terms of its real and 
imaginary parts, xa(w) = imagi¬ 

nary part of the retarded correlation function, Xa('^) i® 
then obtained from the fluctuation-dissipation theorem: 
Xa{^) = (1 ~ at zero-temperature, /? = oo. 


and the second term in parentheses vanishes. The real 
part xSi('^) i® then obtained by the Kramers-Kronig re¬ 
lation: 


xa(^)= r 



TT Uj' — UJ 


(44) 


By choosing A = Zp^a we can measure the local dy¬ 
namical spin density correlations. Explicitly we find that 

{Zp,a{^)Zp, y{t)) = 4(np,<^(0)np/,^/(<)) (45) 

-2np_cr(0) - 2npiy{t) + 1. 

Since npr^a->(t) = np/_o-'(0) the last three terms do not de¬ 
pend on t and only contribute a constant, which does not 
give any contribution to the interesting finite-w structure 
factor. Local dynamical charge and spin correlations are 
straightforwardly calculated from the spin density corre¬ 
lations. 

Similarly, local single-particle excitations can be mea¬ 
sured by choosing 


A = + Cp^a = Zq^aXq^^. (46) 

q<p 


Whether a particle or a hole excitation is realized in 
the excited state that is picked out by phase estimation 
can easily be determined by measuring the total particle 
number of that state. 

Using circuits similar to the preparation of Slater de¬ 
terminants, we can measure dynamical correlation func¬ 
tions not only locally but for arbitrary single-particle 
wave functions. In particular, using momentum eigen¬ 
states one can obtain the momentum-resolved electron 
and hole spectral functions as measured in angle resolved 
photo emission (ARPES). 


VII. NON-DESTRUCTIVE MEASUREMENTS 

For the small systems that we have explicitly simu¬ 
lated, the time needed to adiabatically prepare a good 
estimate of the ground state is small compared to the 
phase estimation time. However, this will change as one 
increases the system size, and the time to adiabatically 
prepare the ground state will ultimately dominate, at 
least if one uses an annealing path that linearly interpo¬ 
lates between initial and final Hamiltonian. The reason 
is the time to adiabatically prepare the ground state is 
expected to scale as for a system with spectral gap 
A, while resolving energies to accuracy A with phase es¬ 
timation takes time 1 /A and thus for small A the phase 
estimation is faster. We note that the improved higher- 
order annealing paths in V C do alleviate this problem. 

This problem is slightly alleviated by the fact that we 
are content to prepare by annealing a state with signif¬ 
icant overlap (say, 1/2 overlap) with the ground state, 
as then phase estimation will give us a significant chance 
of projecting onto the ground state. However, once we 
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have prepared a ground state ipQ by a combination of an¬ 
nealing and phase estimation, it will be desirable to mea¬ 
sure properties of the state without destroying it. Here 
we present two approaches for such non-destructive mea¬ 
surements. Section VIIA discusses an approach based on 
the Hellman-Feynman theorem, which not only is non¬ 
destructive but also scales better in the error than the 
simple way of repeatedly preparing the ground state and 
measuring . An alternate approach based on “recover¬ 
ing” the state after measurements may be more useful in 
some circumstances and is discussed in Sec. VIIB. 


A. Hellman-Feynman based approach 

According to the Heilman Feynman theorem the 
derivative of the ground state energy Eq{X) with respect 
to a perturbation AO is just the expectation value of the 
operator O in the ground state |'I'o(A)): 

J;i?o(A) = ^(d/o(A)|H + AO|vI/o(A)) 

= (vI/o(A)|0|d/o(A)). (47) 

An alternative and superior way of measuring expecta¬ 
tion values of arbitrary observables O is to adiabatically 
add a small perturbation AO to the Hamiltonian H. 

This opens a way to measuring observables without 
destroying the ground state wave function |4>o(0)) of H. 
We adiabatically evolve the ground state wave function 
|4>o(0)) to |4 'o(A)) by slowly increasing A to its final 
(small) value, and then perform another quantum phase 
estimation to determine the energy Eq(X). The ground 
state expectation value can be estimated through 

(4/0(0)1014^0(0)) = J^Ao(A)U.o = ~ +0(A^). 

(48) 

For example, to measure a Green’s function we add a 
“hopping” term A(cj) -|-cJ_£,.Cp_cr) to the Hamiltonian. 

In order to obtain an estimate with an error e we need 
to choose A = 0(e), and then measure the energy to an 
accuracy Ae = O(e^). Phase estimation then requires 
time 0{e~^), which scales the same as repeated prepara¬ 
tion with destructive measurements and thus only gets a 
constant improvement. 

The scaling can be improved by using a symmetric 
estimator for the derivative 

(4/o(0)|0|4/o(0)) = ^o(^) + 0(A3), (49) 

which reduces the approximation error. To get the en¬ 
ergy and thus allows a larger value of A = ©(e^/^), 
which requires phase estimation only to an accuracy of 
C>(e^/^), and a time scaling as In general, 

using a fc-th order estimator for the derivative requires 
0{k) energy measurements, but with an error scaling as 
0{X^), we can choose A = and require time 

Asymptotically for small e we can thus 


estimate the expectation value non-destructively with an 
effort scaling as obtaining a near-quadratic 

speedup. 

In the Hellman-Feynman approach, there may be no 
need to do an adiabatic evolution. In many cases we 
could simply start in the state 'i/’Oi then apply a phase 
estimation to expectation value of -|- AO, and then 
apply another phase estimation to estimate the ground 
state energy of El. The absolute value squared of the 
overlap between the ground state of H and the ground 
state of H -I- AO is equal to 1 — 0(A^/A^), and so for 
A << A the overlap is close to unity. So, with probability 
close to 1, the energy given by phase estimation of i7-|-AO 
will in fact be the desired ground state energy, and with 
again probability close to 1 the second phase estimation 
will return to the ground state of H. 


B. Non-Destructive Measurements With Quadratic 
Speedup 

This annealing procedure becomes less useful for op¬ 
erators O which require a complicated quantum circuit. 
Suppose we wish to measure the expectation value of 
a time dependent correlation such as O = 
where f > 0 and where S^[t) = exp{iHt)S^ exp{—iHt). 
In this case, the quantum circuit to measure O (given 
below) is much more complicated; thus the Hamiltonian 
H + eO becomes significantly more costly to do phase 
estimation on. 

In this section, we explain a technique to measure ar¬ 
bitrary projection operators in a non-destructive fash¬ 
ion. We then further explain a quadratic speedup of 
this approach. Our initial approach will give many in¬ 
dependent binary measurements with outcome probabil¬ 
ity determined by the expection value of the operator. 
This will require a number of samples proportional to 
the square of the inverse error. On the other hand, the 
improved approach will exploit phase estimation to speed 
this up quadratically (up to log factors). 

As we have explained previously, we have quantum cir¬ 
cuits that implement projective measurements of all the 
operators we are interested in, such as number operator, 
hopping operator, and so on. Further, we have a general 
technique for implementing a projective measurement of 
any unitary that can be controlled. 


1. Recovery Map 

Suppose we wish to implement a projective measure¬ 
ment with k possible outcomes, for some given k. We 
describe these outcomes by projectors Qi,...,Qk, with 
Qi — we refer to this measurement as “mea¬ 

suring Q”. The ground state ■i/’o can be written as a 
linear combination of states in the range of each of these 
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projectors: 


^0 = 

i 


(50) 


The transitions are then governed by a Markov chain 
with transition probabilities 

= (1 - |a,n2 = 1 - 2|a,p + \a,\\ (54) 


where 


Qi4^i 


(51) 


and \(j)i\ = 1. Assume for the moment that we can also 
implement a measurement Pq which projects onto the 
ground state i/'o (we describe how to do this to sufficient 
accuracy below). Then, the algorithm for nondestruc¬ 
tive measurement is: hrst, begin in the ground state. 
Then, measure Q. Then, measure Pq. If the result is 
that one indeed is in the ground state, then we have 
restored the ground state; at this point we can either 
re-measure the projector Q to obtain better statistics, 
or make some other measurement of a different observ¬ 
able. If instead the measurement of Pq reveals that we 
are not in the ground state, we re-measure Q and then 
re-measure Pq. We repeat this process of re-measuring 
Q and re-measuring Pq until we are in the ground state. 
The basic idea is similar to that in Ref. 48. 

The rest of this subsubsection is focused on calculat¬ 
ing the probability of returning to the ground state; most 
readers may wish to skip to the next subsubsection, which 
gives the more important quadratic speedup. The main 
reason for introducing the recovery map in this subsub¬ 
section is that it can be used after the quadratic speedup 
(see below). Let us compute the probability of returning 
to the ground state in this process. If the initial mea¬ 
surement of Q gives outcome i then the resulting state 
is (j)i] since, we have a probability \ai\^ 

of returning to the ground state after measuring Pq. We 
now analyze the case that this measurement of Pq reveals 
that we are not in the ground state; then the resulting 
state is equal to 

\ (52) 

Vl - \atV 

= - -(j)i -r'00. 


up to a phase. Then, the probability that the subsequent 
measurement of Q gives outcome j is equal to 






V-o)!' (53) 


and the resulting state is exactly equal to (pj. Repeating 
this, we find that after every measurement of Q with 
outcome i, the probability that the measurement of Pq 
will return to the ground state is |aip, while if we do 
not return to the ground state the probability that the 
next measurement of Q will give outcome j is given by 
Eq. (53). 


— i,j^i — \^i 


(55) 


To^i = |aip, (56) 

where is that, starting in state pi, a measurement 
of Pq reveals that one is not in the ground state, and 
then a measurement of Q gives outcome j, while To<_i is 
the probability that, starting in state pi, does leave on 
in the ground state. The initial measurement of Q at the 
start of the Markov chain gives state pi with probability 
loip. We would like to work out the expected number 
of measurements before the algorithm terminates back in 
the ground state. 

To do this calculation, we re-interpret the probabili¬ 
ties, saying that with probability 2|aip the system makes 
some transition; half the time this transition leads to 
state 0, while the other half the time the system makes 
a transition from state i to some state j, and state j is 
chosen with probability proportional to state \aj\'^. See 
the factor —2|aip in Eq. (54) which we interpret as the 
negative of the probability of there being a transition; 
note also that j may equal i, so some of the “transitions” 
do not actually lead to a change in state: this is the fac¬ 
tor of \ai\^ in Eq. (54). This re-interpretation only makes 
sense if < 1/2, of course but it will allow us to com¬ 
pute the expected time exactly and as the expected time 
is analytic, this calculation will then work for all |aip. 
The advantage of this re-interpretation is that if a tran¬ 
sition occurs, and if the transition leads to a state other 
than 0, then the probability that that state is equal to 
J is proportional to |ajp, which is the same probability 
distribution as the chain began in. 

Starting in state i, the expected number of steps before 
a transition is I/2|aip. Thus, averaging over states i with 
initial probability distribution |aip, the expected number 
of steps before a transition is k/2, where k is the number 
of possible outcomes of measurement Q. Half of these 
transitions return to state 0, while the other half do not, 
leading to some state j 0 with probability distribution 
proportional to |a_, p. Thus, the average number of steps 
before returning to the ground state is 

P^steps — k. (57) 

Note that both the number of phase estimations required 
and the number of measurements of Q required are equal 
to k. 


2. Quadratic Speedup 

Suppose we wish to measure a projector Q. Define 
unitaries U = 2Q — I and V = 2Pq — 1. Assume that Pq 
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has rank one. Applying any sequence of projections Pq 
or Q to the ground state i/jq gives some state in the space 
spanned by the non-orthogonal basis vectors i/'o j QV’o ■ 
Hence, we will work in this two-dimensional space for 
our analysis (the degenerate cases that QV'o = V’o or 
QipQ = 0 can be handled easily too and one can verify 
that the final result will be correct in these cases also). In 
this subspace, in the non-degenerate case, we can write 

2cos 2(6I)-1 -2cos(6»)sin(6l)\ 

y—2cos(0) sin(0) 2sin^(0) —1 y’ 

where the angle 9 is such that (V'olQl^o) = cos^(0). We 
emphasize that the above two equations are written in 
the orthogonal basis ilJo,Z{Qipo — cos^ {0)ipo) , where Z 
is some normalization factor. Thus, the unitary UV re¬ 
stricted to the subspace has eigenvalues 

exp(±2i0). (60) 

We now build a quantum circuit that adds an addi¬ 
tional ancilla that controls both U and V, hence controls 
whether or not UV is applied. We then implement phase 
estimation using this ancilla to determine the eigenvalues 
of UV. We use the ground state ipo as input to the phase 
estimation circuit. Let Tg be the time required to imple¬ 
ment UV. Using the ancilla to control (C/U)", for various 
powers of n, one is able to measure the eigenvalue of UV 
to accuracy e in a time proportional to Tge”^ log(e“^); in 
this case, we take n of order e~^. 

In fact, UV has two distinct eigenvalues, and the phase 
estimation will return one of the two answers randomly. 
However, since the eigenvalues differ only in sign, and the 
desired expectation value {ijjo\Q\fpo) is equal to cos^(0), 
both eigenvalues give the same answer for this expecta¬ 
tion value. 

After implementing this phase estimation to the de¬ 
sired accuracy, we then are left with some state in the 
given subspace. If we then wish to recover the ground 
state (to recover some other observable), we can apply 
the recovery map above by alternating measurements of 
Q and Pg until the ground state is restored. 

3. Measuring Pq 

The above procedure relies on the ability to measure 
Pq. The key is to define an operation that will only iden¬ 
tify whether or not one is in the ground state, without 
revealing additional information. If instead the measure¬ 
ment gives several bits of information on the energy of the 
state, the procedure does not work; intuitively, each (pi 
is a coherent superposition of different energy eigenstates 
and making this more detailed measurement destroys this 
information. 


We now describe how to do this; we call this proce¬ 
dure a “coherent phase estimation”. We have the ability 
to perform a phase estimation on the Hamiltonian H to 
determine the energy of a state. This phase estimation is 
given by a quantum circuit, in which several control bits 
are initialized in a state |0). Then, a Hadamard is ap¬ 
plied to each control bit. Then, the control bits are used 
to control the application of the Hamiltonian for a set 
period of time. Finally, the Hadamard is re-applied, and 
the control bits are measured in the computational basis. 
Then, classical post-processing is applied to extract the 
energy of the state. This process can be performed in se¬ 
rial (using only one control bit) or in parallel; the parallel 
approach uses more control bits but reduces the depth of 
the circuit. The measurement of the energy is not de¬ 
terministic; however, one can employ a larger number of 
control bits to increase the accuracy. 

To measure Pg, we use this phase estimation in its 
parallel form, but do not measure the outcomes of the 
control bits. We initialize an additional “outcome” bit 
to |0). We then use a unitary quantum circuit to im¬ 
plement the classical postprocessing done to determine 
the energy, storing that energy value coherently. Then 
a unitary quantum circuit determines if that energy is 
near the ground state value (which we have determined 
in advance by a phase estimation before doing any mea¬ 
surements), and if so, it flips the value of the outcome 
bit. Finally, we uncompute, reversing the classical post¬ 
processing and phase estimation steps, and then at the 
end we measure the outcome bit. This procedure is es¬ 
sentially that in Ref. 48. While this procedure does not 
implement a projective measurement, but rather imple¬ 
ments a POVM, if sufficient numbers of control bits are 
used, the measurement can be arbitrarily close to the 
desired projective measurement. 

Performing this measurement of Pg does require addi¬ 
tional qubits. The number of qubits required scales pro¬ 
portional to the number of bits of accuracy desired; how¬ 
ever, since this accuracy is of order A, if A is only poly- 
nomially small in system size then this requires only log¬ 
arithmically many extra bits. Additionally, the more ac¬ 
curately we implement a projective measurement rather 
than a POVM, the more control bits required, but this 
overhead also only scales logarithmically. 

An alternate method of implementing Pg is to note that 
we can easily implement the projector onto the ground 
state of the Hamiltonian at the start of the annealing pro¬ 
cess if our initial Hamiltonian is of a sufficiently simple 
form, such as a free fermion Hamiltonian or a Hamil¬ 
tonian with a product ground state. For a general free 
fermion Hamiltonian, we can use Givens rotations to ro¬ 
tate to the case that the ground state simply has some 
spin-orbitals occupied and some empty; we then mea¬ 
sure whether we have the desired pattern or not. Call 
this projector Pq^ . Then, if Uadiab is the unitary which 
describes the approximately adiabatic evolution from the 
free Fermion Hamiltonian onto the desired final Hamil¬ 
tonian, we can approximate Pg « UadiabPl^ The 
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accuracy of this approximation depends upon the anneal¬ 
ing path; see the discussion in V C, which shows that 
we can use this strategy to achieve a (near) quadratic 
speedup. The improved annealing paths discussed there 
are most useful for the particular approach here where we 
need very high accuracy in the annealing; in all other ap¬ 
plications in the paper, we need only moderate accuracy 
in the annealing, sufficient to get a large overlap with the 
ground state so that phase estimation can then project 
onto the ground state with reasonably large probability. 
We therefore 


VIII. CONCLUSIONS 

In this paper we have presented a comprehensive strat¬ 
egy for using quantum computers to solve models of 
strongly correlated electrons, using the Hubbard model 
as a prototypical example. Similar strategies can be used 
for generalizations of the Hubbard model and other dis¬ 
crete quantum lattice models, including, but not limited 
to, the t-J model, frustrated magnets, or bosonic mod¬ 
els coupled to (static) gauge helds, which all suffer from 
negative sign problems in quantum Monte Carlo simu¬ 
lations on classical computers. We go beyond previous 
papers'^’^^^^® that discussed simulations of the Hubbard 
model on quantum computers by giving complete details 
of all steps needed to learn about the properties of the 
Hubbard model. 

In particular we presented a strategy to prepare trial 
ground states starting from various mean-held states 
or RVB states on decoupled plaquettes. While hnding 
the ground state of quantum lattice models is QMA- 
hard'^’'^®’^° this is a statement about the worst case. We 
expect that, as experience has shown, many experimen¬ 
tally relevant models have ground states with either no 
long-range order (so-called quantum liquids) or with bro¬ 
ken symmetries that are qualitatively well described by 
mean-held theories. In fact, if a model has such pecu¬ 
liar properties that it is hard to hnd its ground state 
on a quantum computer, then also a material described 
by that model will have a hard time reaching its ground 
state. Quick low-temperature thermalization in experi¬ 
ments is thus an indication that the ground state of the 
model describing its properties should also be easy to pre¬ 
pare on a quantum computer. The opposite is the case, 
e.g. for (classical) spin glasses, where hnding the ground 
state is NP-hard but also spin glass materials never ther- 
malize nor reach the ground state. 

We have furthermore presented an efficient and deter¬ 
ministic quantum algorithm to prepare arbitrary Slater 
determinants as initial state, which scales better than the 
algorithms of Refs. 5 and 15. This can be used to prepare 
ground states of various candidate mean-held Hamiltoni¬ 
ans, from which an adiabatic evolution to the ground 
state of the Hubbard model can be attempted. 

To implement time evolution under the Hubbard and 


mean-held Hamiltonians we have given explicit quantum 
circuits for all terms, and discussed the size of the Trot¬ 
ter time step required to achieve sufficiently small er¬ 
rors. We discuss how time evolution under the individ¬ 
ual terms in the Hubbard model can be efficiently paral¬ 
lelized, ultimately requiring 0{N) qubits and gates with 
only O(logV) parallel circuit depth for one time step, 
which allows efficient simulations of very large systems. 
We gain additional constant factors over previous ap¬ 
proaches, by optimizing the phase estimation algorithm 
to reduce the required number of (expensive) rotation 
gates by a factor of four, assuming that our gate set 
consists of one- or two-qubit Clifford operations and ar¬ 
bitrary single qubit rotations. We furthermore propose 
to use a larger Trotter time step for the adiabatic state 
preparation (whose time scale is controlled by the inverse 
gap squared), to prepare a very good (but not perfect) 
guess for the ground state, and then rehne this to the ex¬ 
act ground state by doing only the final quantum phase 
estimation (whose time scale is shorter since it is con¬ 
trolled by the inverse gap) with a small time step. 

We have finally discussed approaches to obtain a 
quadratic speedup in measurements by proposing two 
non-destructive measurement strategies, one based on 
the Hellman-Feynman theorem, and another based on 
recovering the ground state after a (destructive) mea¬ 
surement of a single qubit. We also introduce a new 
approach of measuring dynamic structure factors and 
spectral functions directly in frequency space, using ideas 
similar to angle-resolved photoemission experiments. 

Estimating the gate counts required to simulate the 
Hubbard model, we find that even on lattices with more 
than N « 1000 sites - which should be large enough to 
learn most interesting properties - can be simulated on 
small scale quantum computes using only a few thousand 
qubits and a parallel circuit depth of about one million 
gates. This number is based on our estimate of a few 
thousand time steps with a circuit depth of not more 
than a few hundred gates each. This means that even 
with logical gate times of the order of l^s, the ground 
state of the Hubbard model can be prepared within a few 
seconds. Small scale quantum computers will thus be a 
very powerful tool for the investigation of many problems 
in the held of strongly correlated electron models that are 
currently out of reach of any classical algorithm. 
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Appendix A: Error bounds for time-dependent 
Trotter Formulas 

Simulating the dynamics of time-dependent Hamil¬ 
tonians on quantum computers is a more subtle issue 
than simulating time-independent dynamics and results 
proven for the time-independent case do not necessar¬ 
ily transfer over. This issue is significant here because 
of issues surrounding simulating adiabatic state prepara¬ 
tion. Bounds for the error in high-order Trotter formu¬ 
las are known, however, the scaling of these bounds with 
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the number of terms in the Hamiltonian is known to be 
loose. Since the Hamiltonians needed for our purposes 
have a large number of terms, tighter estimates of the 
error scaling may be important for determining the cost 
of performing adiabatic state preparation for Fermionic 
systems. 

It may be tempting to think that we can neglect the 
errors in the Trotter formula that are incurred by ap¬ 
proximating H{t) by a piecewise constant Hamiltonian 
because the time derivatives of the Hamiltonian are small 
for a slow adiabatic evolution, but this is not necessar¬ 
ily true because adiabatic theorems only require that the 
derivatives of the Hamiltonian are small relative to an 
appropriate power of the minimum eigenvalue gap. This 
means that in some circumstances adiabatic evolution 
may actually involve a rapid passage. Furthermore, since 
we are interested in high accuracy state preparation even 
though such errors may be small they will not necessarily 
be negligible compared to our target error tolerance. 

In this section we analyze these errors in more detail. 
In particular, we will examine the error in the second 
order order Trotter-Suzuki formula 


m 1 

i=l j=m 

where we approximate the evolution from time 0 to 1 by 
a single second-order Trotter-Suzuki step. The bounds 
that we present are a generalization of those in Ref. 10 to 
the time-dependent case. Bounds for the Trotter formula 
are given in Ref 51. 

Our main result is that 
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This result reduces to the bound of 10 and agrees with 
the asymptotic error scaling predicted in 52, up to a 
use of the triangle inequality and a small multiplicative 
factor, if H is time-independent. 

We prove the bound in two steps. First, we will show 
that 
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Then, we use the bound of Ref. 10 that 
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Eq. (A2) then follows from the triangle inequality. We 
remark that in Ref. 10, a typographical error led to 
\\[[A,H{x)],H{x)]\\ < ||[[A,i?],A]|| + ||P,H],H]|| be¬ 
ing replaced by \\[[A,H{x)],H{x)]\\ < ||[[A,H],A]] -b 
[[A,i3],H]|| in the language of that paper in appendix 
B (see text above Eq. (B3) of that paper). This pair of 
missing || symbols propagated to later works although it 
does not affect any of the scaling results cited in later 
papers. In the above equation, we have corrected this 
error, so that the right hand-side is a sum of two distinct 
norms for each j, rather than a norm of a sum. 

To show Eq. (A3), we write to denote the 

first, second, ... derivatives of H with respect to time. By 
Taylor’s theorem, H{u) = H{l/2.) + {u — l/2)H'{l/2) + 
Ii/ 2 ^'^ — s)Fl"(s)ds. Hence Te~^-l^o ^(“)d« jg 


Te 
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We next bound the difference 
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Transforming to the interaction representation gives 

(^H(l/ 2 ) + (u-l/ 2 )H'(l/ 2 )'^ du 
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where we define 
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So by a unitary rotation 
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Note that ||i7^(l/2) — 77'(1/2)|| can be written as 






which can be used to show that ||i7^(l/2) — iJ'(l/2)|| < 


|u-l/2|||[iJ'(l/2),77(l/2)]||,and so 

< f \u-ll2\^[H'{l/2),H{l/2)]\\du 
Jo 

= l||[i7'(l/2),i7(l/2)]||. (A7) 

Eq. (A3) then follows from Eqs. (A5,A6,A7) and the 
triangle inequality. 






